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NUMERICAL  COMPUTATION  OF  HYPERSONIC  FLOW  PAST 
A  TWO-DIMENSIONAL  BLUNT  BODY 


ABSTRACT 
Given  a  two-dimensional  blunt  body  in  a  steady  super- 
sonic flow  of  compressible  fluid,  the  problem  is  to  study  the 
flow  between  tne  body  and  the  detached  shock  wave  formed  in 
front  of  it.   In  this  paper,  the  shock  shape  is  assumed  to  be 
analytic  and  we  use  it  as  an  initial  line  in  solving  a  Cauchy 
problem  for  tne  flow  and  the  body.   Hyperboliclty  of  the  partial 
differential  equations  of  motion  Is  achieved  in  the  subsonic 
region  by  means  of  complex  extension  on  one  of  the  Independent 
variables.   By  means  of  such  a  substitution  the  first  order 
system  of  nonlinear  equations  in  two  dimensions  is  transformed 
into  a  symmetric  hyperbolic  system  In  three  independent  vari- 
ables.   Such  an  Increase   in  the  number  of  independent  vari- 
ables can  be  avoided  if  the  equations  are  written  in  normal  form 
with  respect  to  two  characteristic  coordinates ,   Difference 
Schemes  of  second  order  accuracy  are  derived  for  both  formula- 
tions of  the  complex  extension  and  chey  are  programmed  for  the 
IBM  709^0   Five  digit  accuracy  is  obtained  In  the  solution  and 
numerical  results  are  presented  for  a  particular  flow  at  Mach  6, 
The  relative  merits  of  the  two  formulations  are  discussed.   Use 
of  the  three-dimensional  symmetric  hyper'bolic  system  is  simple 
because  the  equations  are  valid  In  the  entire  domain  of  solution, 
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However,  stability  considerations  restrict  the  mesh  ratios  so 
that  the  body  Is  barely  reached.   Less  trouble  Is  encountered 
In  solving  past  the  body  when  the  equations  In  normal  form 
are  used.   A  serious  drawback  In  this  second  approach  is  that 
the  characteristics  cusp  at  the  sonic  line.   Thus,  not  only 
is  one  prevented  from  integrating  across  the  sonic  line,  but 
also  a  portion  of  the  flow  is  sometimes  made  inaccessible  when 
characteristic  coordinates  are  used. 
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NUMERICAL  COMPUTATION  OF  HYPERSONIC  PLOW  PAST 
A  TWO-DIMENSIONAL  BLUNT  BODY 

by 

Eva  V.  Swenson 
1.    Introduction 


In  this  paper  we  shall  consider  a  two-dimensional 
blunt  body  in  a  steady  supersonic  flow  of  compressible  fluid. 
Our  aim  is  to  study  the  behavior  of  the  flow  between  the  body 
and  the  shock  wave  caused  by  the  body.   The  direct  approach 
consists  in  prescribing  the  body  shape  and  solving  the  partial 
differential  equations  of  motion  of  the  flow.   However,  the 


shock- 


streamlines 


subsonic  region 


sonic  line 


supersonic  region 
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problem  Is  complicated  by  the  fact  that  the  location  and  shape  of 
the  shock  wave  determined  by  the  given  body  are  unknown.   An  alter- 
native approach,  referred  to  as  the  Inverse  method,  is  to  start 
from  a  known  shock  shape  and  to  find  the  flow  and  the  generating 
body  by  solving  a  Cauchy  problem  for  the  equations  of  motion.  In 
the  region  where  the  flow  Is  supersonic,  the  Cauchy  problem  is 
properly  set  because  the  equations  of  motion  are  of  hyperbolic 
type.   However,  near  the  vertex  of  the  shock  wave,  where  it  is 
known  that  the  flow  is  subsonic,  the  partial  differential  equations 
of  motion  are  elliptic  and  initial  data  on  the  shock  wave  are  not 
appropriate  to  determ.lne  a  solution. 

Various  methods  have  been  applied  to  overcome  the  diffi- 
culty just  described,  including  the  use  of  double  power  series 
expansions  and  smoothing  processes;  see  Lewis  [10],  and  Hayes 
and  Probstein  [8],  which  contain  a  good  exposition  of  the  work 
to  be  done  in  these  directions.   In  our  case,  we  choose  Instead 
to  follow  the  procedure  described  by  Garabedlan  [6,7],  which  ex- 
ploits the  analyticity  of  the  problem  to  continue  the  equations 
into  the  complex  domain.   The  use  of  complex  extension  on  one  of 
the  independent  variables  leads  to  hyperbolicity  of  the  equations 
and  the  Cauchy  problem  becomes  properly  posed.   One  version  of 
this  method  based  on  characteristic  coordinates  has  been  applied 
to  the  problem  of  three-dimensional  flow  past  a  blunt  body  of 
revolution  [7].   For  that  case,  however,  difficulties  were  experi- 
enced along  the  sonic  line. 

For  two-dimensional  flow,  we  shall  present  in  Sections 
2  and  3  three  possible  formulations  of  the  problem,  the  third  of 


-  6  - 


which  Is  a  first  order  system  with  the  velocity  components   u,  v 
and  the  stream  function  i{/      as  dependent  variables.   We  find  that 
the  latter  formulation  eliminates  the  difficulties  previously  en- 
countered in  solving  for  the  flow  in  the  transonic  region.   In 
Sections  4  and  5,  complex  extension  is  applied  to  the   u,  v,  ip- 
version  in  its  two  forms  -  both  as  a  first  order  system  and  as  a 
system  in  characteristic  normal  form. 

Difference  schemes  are  devised  for  the  numerical  solution 
of  the  two  formulations  of  the  problem,  and  their  stability  is 
discussed  in  Section  6.   Both  of  these  schemes  have  been  programmed 
for  the  IBM  7O94,  and  in  Section  7  we  evaluate  the  numerical  re- 
sults for  a  specific  flow  at  Mach  6.  The  stability  and  convergence 
of  the  numerical  scheme  in  practice  is  demonstrated  in  the  tables. 
We  have  found  that  an  accuracy  of  five  significant  figures  can  be 
achieved  readily  in  the  numerical  calculations.   Moreover,  the 
results  are  in  satisfactory  agreement  with  the  experimental  data 
presently  available  on  two  dimensional  flows  past  blunt  bodies, 
see  [1,  9,  11]. 

The  author  would  like  to  express  deepest  gratitude  to 
Professor  Paul  R,  Garabedian  who  suggested  the  problem  and  who  has 
been  an  invaluable  source  of  inspiration,  encouragement,  and  advice. 

2.    Equations  of  Motion 

We  shall  start  with  the  equations  of  motion  for  steady 
plane  flowT  of  a  perfect  compressible,  inviscid  fluid  as  derived 
in  [5].   Let   x,y  be  Cartesian  coordinates  in  the  plane.   We 
shall  denote  the  horizontal  and  vertical  components  of  the  velocity 


-  7  - 


by  u(x,y)   and  v(x,y)  ,  the  pressure  by  p(x,y)  ,  and  the 
density  by   p(x,y)  . 

That  mass  is  conserved  is  expressed  by  the  continuity 
equation 


(2.1)  ^4^  +  ^4^^  =  0 

ox      oy 


Since  pressure  is  the  only  force  acting  on  the  fluid,  the  con- 
servation of  momentum  in  the   x   and  y   directions,  respectively, 
is  expressed  by 

(2.2)  u^  +  v|H+I|£=0, 

dx  oy        p      dx 

dx  dy        p      dy  ' 

From  the  law  of  conservation  of  energy  we  get  the  fourth  equa- 
tion of  motion 


2.4)         .u  1^  +  ^  |S 

dx     dy 


which  states  that  if  the  work  done  is  due  solely  to  the  pressure 
forces,  then  the  specific  entropy   S   at  a  moving  particle  remains 
constant.   If  we  assume  that  our  fluid  satisfies  the  equation  of 
state  for  a  polytroplc  gas 

(2.5)  P  =  A(S)p^   ,        7  constant  , 


which  asserts  that   pp  '^   is  a  function  of  entropy  alone,  then 
we  can  rewrite  (2.'^)  as 


2.6)  u  1^  (pp-^)  +  V  |p  (pp-^)  =  0 


The  constant  j      is  the  ratio  of  specific  heats  of  our  fluid. 

Equations  (2.1),  (2.2),  (2-3)  and  (2.6)  together  with 
relevant  boundary  conditions  are  sufficient  to  define  the  motion 
of  a  compressible  fluid.   We  shall  often  use  the  conservation 
of  energy  equation  in  the  integrated  form 


(2.7)  \   (u^  +  v^)  +  -L—     ^   =  H   , 


where   H  is  a  constant.   Equation  (2.7)  is  referred  to  as 
Bernoulli's  equation, 

A  more  convenient  form  of  the  equation  of  state  (2.5' 
can  be  arrived  at  if  we  note  that  the  continuity  equation  (2.1' 
implies  the  existence  of  a  stream  function  ip{y:,y)      such  that 


(2.8)  ^   =  up  ,       -  "^„  =  vp  , 

y  ^ 


that  is,  a  function  ip     which  is  constant  along  each  particle 
path,  or  streamline.   Thus,  we  can  rewrite  (2.5)  as 


(2.9)  P  =  A(^)p^  , 


where   A   can  be  considered  as  an  arbitrary  function  of  f      alone, 


Our  next  task  Is  to  pick  a  suitable  set  of  equations 
for  our  problem.   The  most  obvious  choice  Is  the  first  order 
system  for   u,v,p,p   consisting  of  (2.1),  (2.2),  (2.3)  and  (2.6), 
which,  written  In  matrix  form.  Is 


(2.10) 


^ul 

'  uvp 

D 

7P 
D 

V 

D 

0 

u 

V 

0 

V 

u 

1 

up 

0 

V 

p 

7vpp 

D 

7upp 
D 

uvp 

D 

0 

p 

p 
V  ) 

X 

D 

1^ 

2 

uD 

V 

"  7 

p 

y 


where 


2       2 
D=7p-up=p(c 


u 


and   c  =  ydp/dp   Is  defined  as  the  speed  of  sound.   There  Is 
a  considerable  disadvantage  In  using  this  system.   Note  that  It 
does  not  take  Into  account  the  fact  that  the  conservation  of  energy 
equation  (2.6)  can  be  Integrated  directly  to  give  Bernoulli's 
equation  (2.7).   Not  only  does  (2.6)  add  to  the  order  of  the 
system  (2.10)  unnecessarily,  but  It  also  Introduces  the  stream- 
lines as  a  family  of  characteristics.   This  is  not  desirable 
because  near  the  body  the  streamlines  bend  in  a  direction  parallel 
to  the  body  and  perpendicular  to  the  direction  in  which  we  wish 
to  Integrate  the  system  (2.10).   We  can  never  expect  to  get  near 
the  body  with  this  formulation.   Hence  we  intend  to  retain  this 
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version  only  for  checking  purposes « 

A  second  way  of  formulating  the  problem  is  motivated 
by  the  desire  to  solve  for  f      expllclty.   The  stream  function 
Tp      Is  of  Interest  to  us  because  the  streamline  f   =  0     will 
serve  to  locate  the  body  generating  a  given  shock  wave.   From 
equations  (2,2),  (2.3),  (2-7),  and  (2-9)  one  can  derive  (see 
[3,7])  the  following  second  order  quaslllnear  partial  differ- 
ential equation  for  ^  : 


^om^   /2   2,,    ^    ,   ^,2  2.  ,       ^   2    2    k'  [Tp]       2H+U  +v     ^ 

2.11     c  -u  )^   -2uv  i^      +(c  -V  )^   +p  c    .  )  ,    or =  0 

'  ^xx     ^xy         yy       A(^^)      27 


Together  with  equations  (2.8)  and  Bernoulli's  equation  (2.7) 
In  the  form 


(2.12)     I  p-2  iijl    +   i'p    +  ^      M-P)    p^"^  =  H   , 


(2.11)   gives  an  alternate  equation  for  our  problem.   From 
equation  (2.11)  and  (2.8),  a  first  order  system  of  partial 
differential  equations  In  the  dependent  variables  ip   ,    f        and 
Tp        can  be  obtained. 

y 

Superficially,  It  might  seem  as  though  we  have  re- 
duced our  problem  to  that  of  solving  three  differential  equa- 
tions instead  of  four.   However,  the  work  of  integrating  a  fourth 
differential  equation  has  been  replaced  by  that  of  solving  for 
p   as  a  root  of  Bernoulli's  equation  (2.12).   By  applying  Newton's 
method,  we  can  obtain  p   to  second  order  accuracy.   However, 
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2     2     2 
at  the  sonic  line,  we  have   u  +  v   =  c    and  therefore 


^p  =  "P"^  ^^x  ^  '4^  ^  7Ap^"2 


2,2     9 
^   +  ^   +  ^  =  0 


so  that  Newton's  method  is  no  longer  applicable  and  we  have  to 

resort  to  more  sophisticated  schemes. 

In  order  to  circumvent  this  difficulty,  we  introduce 

the  set  of  dependent  variables   ^  ,  u  ,  v   instead  of  f    ,    f      ,    ip 

^    y 

This  brings  Bernoulli's  equation  (2.12)  into  a  form  that  is 
solvable  for   p   explicitly  as 


P  = 


u  +  V   V   7-1  -ji-r 


(H  -  "  :  ^   )  -.r^^iv] 


2     '      7k[ilj) 


We  shall  proceed  to  derive  a  set  of  equations  for   u  ,  v  ,  and  ifj 
based  on  the  latter  relation. 

If  we  differentiate  (2-7)  with  respect  to  x   and  y  , 
we  get  the  equations 


2    A!  2 

(2.13)       uu^  +  vv^  +  -^  ^  lb      +  —  p   =:  0   , 
X     X    7-1   A   ^x    p   ^x       ' 


2  2 

(2.1^)       uu   +  vv^^  +  -^  ^  ip      +  —  P   =  0 
y     y   7-1  h     ^Y       P   y 


Multiplying    (2.13)    by     u    ,    (2.14)    by      v    ,    and   adding,    we    obtain 


(2.15)  u   u      +   uvv      +  uvu      +  v^v      +  —   (up      +  vp    )    =  0      , 

X  X  y  y        p  ^x  ^y'  ' 
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where  (2.8)  determines  Tp        and  f      .      Using  (2.1)  we  then  re- 

X        y 

write  (2.15)  as 


(2.16)       (c^  -  u^)  u^  -  uv(u^  +  v^)  +  (c^-  v^)v^^  =  0 


This  Is  the  first  equation  of  our  system. 

The  second  equation  may  be  obtained  by  subtracting 
(2.2)  from  (2.13)  to  get 


?  2 

p  A  '  f         1 

2.17)       vv   -  vu   +  -^-  ^  Tp      +-^^p P^  =  0   . 

'  '        X     y   7-I  A   ^x  p   "^x    p   X 


We  can  simplify  (2.17)  if  we  note  that 


P^  =  A'  Tp^p'^    +  yAp^p'^'         ,  c   =  7Ap^ 


Then  (2.17)  becomes 


2  2 

V   -  u    -  — T-  -A~  vp  +  —  "ft—   vp  =  0   , 
X     y  y  -1      A         '^         y         A  ^  ' 


or 

(2-18)     v^  -  Uy  =y|frry  I^P 


We  would  like  to  remark  that  one  can  derive  this  same  equation 
by  working  with  (2.14),  (2.3)  and 

Py  =  A'^yP^  +  rApyP^"^ 


-^  - 


Instead.   To  complete  the  set,  we  may  use  either  of  the  equations 


(2.8: 


f^    =   -vp  , 


^y   =   l^P      > 


depending  only  on  our  choice  of  the  distinguished  Independent 
variable.   Thus  altogether,  our  system  In  the  unknowns   u,  v,  ^   Is 

uv 


(2.19) 


'u 


^^ 


/  2uv 

/   2   2 
/   c  -u 


2   2 

C  -V 


0 


2   2 
c  -u 


0 
0 


+ 


2   2 

c  -u 


Q 


■vp 


where 


c^p    A'(^) 
7(7-1)  ~KW 


Taking  Bernoulli's  equation  (2.7)  Into  account,  we  have  a  suit- 
able formulation  of  our  flow  problem. 

We  shall  be  Interested  ultimately  In  a  Cauchy  problem 
for  the  first  order  system  (2.19)  with  data  prescribed  behind  a 
shock  wave  whose  equation  has  the  form  x  =  f(y).   Therefore  let 
us  consider  the  change  of  coordinates 


y 


C  =  X  -  f(y) 


which  has  the  effect  of  mapping  the  initial  curve  onto  the  line 

^  =  0   In  the  t,,    ^-plane.   The  system  of  equations  (2.19)  which 

in  matrix  notation  has  the  form  U   =  RU  +  S   is  then  transformed 

X     y 
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into 


U   =  [I  +  f'(|)R]~^  R  U^  +  [I  +  f'(e)R]  ^  S   , 


or  more  specifically. 


2.20) 


where 


u^ 


Ui 


D 


/2uv  -I  f '  (^)  (c^-v^) 


2   2 
c  -u 


0 


+ 


D 


I    2   2 
-  c  -V 


f'(5)(c^-v^)   0 


0 


Q[uv  +  f'(^)  (c^-v^) 


Q[c^-u^  +  uv  f '(1)] 


-  vpD 


D  =  (c^-u^)  i-  2uv  f'(^)  +  f'^(?J  (c^-v^)   , 


c^    A  '  ( ?// ) 
7(7-1)  "TH^T  P   ' 


This  is  one  of  the  principal  formulations  of  the  equations  of 
motion  that  we  shall  use  for  the  numerical  solution  of  the 
detached  shock  problem. 

A  more  subtle  approach  which  we  shall  have  occasion  to 
exploit  is  the  reduction  of  equations  (2.l6),  (2.l8)  and  (2.8) 
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to  normal  form  (see  [4,6]),  that  is,  a  form  where  each  equation 
Involves  a  differentiation  In  only  one  characteristic  direction. 
We  shall  confine  our  attention  to  equations  (2.l6)  and  (2.18), 
since  (2.8)  can  be  put  in  characteristic  normal  form  easily. 
In  matrix  notation,  (2.l6)  and  (2.l8)  can  be  written  as 


(2.21) 


/c  -u     -uv 


V°      V 


AA 


v 


+ 


/        2   2 
/-uv     C  -V 


\-l 


0 


u 


0 


A  W 


=  0  , 


where 


Q  = 


7(7-1)    A  (7//)  P 


Given  initial  values  of  the  vector  col(u,v)  =  w  on  an 
initial  curve   *(x,y)  =  0  with  {<t>    )      +    {<t>    )      /^  0  ,  we  are  to 
determine  the  vector  w  whose  first  derivatives  satisfy  a  system 
of  the  form  Ew  +  Fw  .  +  J  =  0  .   It  has  been  proven  in  [4]  that 

^    y 

a  necessary  and  sufficient  condition  for  the  unique  determination 
of  all  the  first  derivatives  along  <i>{y:,y)    =0   is  the  non-vanish- 
ing of  the  characteristic  determinant  of  the  system 


(2.22: 


Et  -  F 


2   2 
(c  -u  )t  +  uv 


-  UVT 


2   2. 
c  -V  ) 


T 


^  0  , 


where   t  =  -  't>  /<i>      .      If   T(x,y)   is  a  real  solution  of 
-^  y 

|Et  -  F|  =  0  ,  then  the  curves  along  which   t  -=   dy/dx  =  const 
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are  called  characteristic  curves.   In  our  case,  the  characteristic 
slopes  are  solutions  of 


2.23)       (c2-u2)  (§|)2  +  2UV  ^+  (c^.v^)  .  0   , 


that    is. 

dx    " 

^   ^+   - 

(2.24) 

/  2      2      2 
-UV+    c  v/u    +v    -c 

2        2 
C     -U 

If  we  solve  the  system  of  linear  homogeneous  equations 

A(Et  -  F)  =0    for  the  vector   A  =  (A-,,A  )  ,  our  original  system 

Ew   +  Fw   +  J  =  0   can  be  rewritten  as  the  linear  combination 
X     y 


AEw   +  AEtw   +  AJ  =  0   , 
X       y 
or  ■^ 


(2.25)  AE(w   +  TW  )  +  AJ  -  0   , 

X     y  ' 


where  now  the  unknown  w  is  differentiated  along  one  character- 
istic direction. 

Carrying  out  the  above  procedure  in  more  detail,  we 
find  that  the  components   A, ,  Ap   of  the  vector  A   satisfy  the 
ratios 


(2,26)       ^  =  -  [(c2-u2).  +  uv]  =Hvl_±Jci^    ^ 

2   2 
Setting   -^1=1  '"-^^^      ^2  =  ~  t  ( '^  "^  )'^  +  ^"v  ]  ,  we  have  that 
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1\        times  (2.16)  plus   A^   times  (2.l8)  gives 


(2.27) 


c2-u^)u^  +  [(c^-u^)T]Uy  -  [(c^-u^)t  +  2uv]v^  +  (c^-v^)Vy 


=  -  Q[(c^-u^)t  +  uv] 


Using  (2.26),  equation  (2.27)  becomes 


(2.28)   t(c^-u^)(u^+tu  )  +  (c^-v^)(v^+TV J  =  -  Qt[(c^-u^)t  +  uv]  , 


where   r   is  defined  in  (2.24). 

Let   a(x,y)  =  0  and   p(x,y)  =0  be  the  characteristic 
curves  with  slopes   t   and   t  ,    respectively.   Then 

da   ox  a   oy  '^  a    a  ox   x   dy     a  dx    +  oy 


Similarly 


Hence  (2.28)  can  be  rewritten  as 


o      o  p   p  /  p   p   2 

T^(c  -u  )u^  +  (c  -V  )v^  =  -  Q  c  /u  +v  -c    y^ 


(2.29) 


T_(c  -U  )Ug  +  (c  -V  )Vo  =       Q  c  /u' 


2,2   2 


^O'V^-V/Vg    -     ■qiC  /u  +v   -c      Yq 


Equations  (2.29)  together  with 


-  1^ 


(2.30) 


y       -    T      X       =   0    , 


yp  -  ^-  ^p  =  °  ' 


and  (2.8)  written  as 


(2.31)  Vp  +  pv  Xp  -  pu  y^  =  0  , 


comprise  a  system  in  normal  form  equivalent  to  (2.I9). 

3.    Shock  Conditions  and  Initial  Conditions 

The  Inverse  detached  shock  problem  consists  of  starting 
from  a  given  shock  and  solving  the  equations  of  motion  for  the 
flow  and  the  generating  body.   In  other  words,  we  have  to  solve 
a  Cauchy  problem  for  the  system  (2.I9),  or  equlvalently  (2.29  -  31 ) j 
with  Initial  data  on  the  shock  wave.   In  this  section  we  shall 
Indicate  how  initial  data  along  a  shock  wave  are  calculated  from 
the  shock  conditions. 

Across  the  shock,  pressure  and  density  are  discontinuous, 
but  the  conservation  laws  enable  us  to  relate  the  values  of  these 
quantities  on  either  side  of  the  shock.   Adopting  the  notation  of 
[3],  which  contains  a  derivation  of  the  shock  conditions,  we  let 
N  and  L   denote  the  normal  and  tangential  components  of  the  vel- 
ocity vector  with  respect  to  the  shock.   We  shall  use  the  sub- 
cript   "o"   to  indicate  quantities  In  front  of  the  shock,  and   "l" 
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to  denote  quantities  behind  the  shock. 

That  mass  Is  conserved  across  the  shock  Is  expressed  by 

(3.1)  P^N^  =  P^N^  =  m  , 

where   m  Is  the  mass  flux  across  the  shock  surface. 
Conservation  of  momentum  gives  the  equations 


(3.2)  p^n2  +p^  =  p^n2  +p^  =  P  ,    L^  =  L^ 


where   P   Is  the  total  flux  of  momentum  normal  to  the  shock. 
Finally  we  have 


■^o  ^1 


which  states  that  energy  is  conserved  across  the  shock. 

Given  the  quantities   L   ,  N   ,  p   ,  p    in  front  of  the 
^  o    o  '  ^o   "^o 

shock,  we  shall  compute  the  quantities   N,  ,  p,  ,  p.,   behind  the 
shock  using  the  shock  conditions  (3.1  -  3).   Since  these  equations 
are  nonlinear  in  N   ,  p-,  ,  and  p,  ,  we  have  to  use  the  auxili- 
ary equations 


Pn   -  P 

(3-4)  N  N,  =   1     ° 


o  1    Pi  -  Po 


> 


and 


(3-5)  Pi  =  PqN^  (1-^x2)  .  ^% 
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The  first  is  referred  to  as  Prandtl's  relation;  the  second  as 

Hugonlot's  relation.    Both  may  be  derived  by  performing 

a  few  algebraic  manipulations  on  the  shock  conditions  (J.l  -  3)- 

Hence,  given   p   ,  p   ,  N   ,  and   L    In  front  of  the 
J  to       ^o     o     o  o 

shock,  we  can  solve  for   L,  ,  p-,  ,  N-,  ,  and   p-.   In  the  sequence 


L-,  =  L  , 
1    o  ' 


P,  =  Po^od-^") 


2 

M-  P, 


(3.6) 


N^  =  [^^N^  +  (1-^1^)7  ]/Nq   , 


p,  =  p  N  /N^   . 
^1    '^o  o'  1 

In  our  particular  problem  we  shall  take  the  direction  of 

Incident  flow  to  be  parallel  to  the  axis  of  the  body,  which  is 

Itself  parallel  to  the  x-axls.   Thus   v   =0  and  u   =  U   ,  a 
^  o  o    o 

constant  In  front  of  the  shock.   Let  6      be  the  angle  which  the 
normal  to  the  shock  wave  makes  with  the  horizontal.   Then  with 
respect  to  the  shock  wave  x  =  f (y)  ,  we  have 
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(3.7) 


L   =  u   sin  6 
o  o 


N   =  u   cos  0 
o    o 


u^f  (y) 

/i+f'^(y) 

u 

o 

/i+f'^(y) 


without  loss  of  generality,  we  shall  choose  a  scale  so  that 

p   =  p   =  1  .   The  speed  of  a  flow  is  usually  specified  by 
■^o    o 

giving  its  Mach  number  M  defined  as 


(3.8) 


M  = 


/  2^  2 
,/u  +v 


where  ^u'^+v'^      is  the  speed  of  the  flow  and   c   is  the  local 


2,  2 


speed  of  sound.   In  our  case   c   =  7P  /p  =  7    ,    hence 

o      o  "^o 

u   =U  =   c    n   =    Jy    W    . 
o  o  o  ^  ' 

Using  the  equations  (3.6)  we  then  have  expressions  for 
L   ,  N,  ,  p   ,  p.,   compatible  with  our  conditions  on  the  incident 
flow.   They  are  as  follows: 


L. 


/y    M   f'(y) 

v/i+f'^(y) 


P- 


(l-^x^)    y    M^  2 

i+f'^(y) 


(3.9: 


N. 


fgy  +  tiLxJlLl       /i+f'^(y)     y 
1'^+^       l+f'^(y)J  /7     M 


7M 


2 


i+f'^(y) 


2y        M.^   y    M^ 

ly+i       .  ,.,2, 


i+f '    (y)J 
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Finally,  we  obtain  u.,  ,  v   from  the  relations 


(3.10) 


u,  =  N-,  cos  e  +  L-,  sin  9    , 


V   =  L   cos  0  -  N   sin 


At  this  point,  we  would  like  to  note  that  the  expres- 
sions for   U-,  ,  v^  ,  p,  and   p-,   consist  of  elementary  and  an- 
alytic formulas  depending  only  on   f'(y).   They  are  elementary 
In  the  sense  that  only  the  operations  of  addition,  multipli- 
cation, division  and  square  root  are  Involved.   If  the  expres- 
sion f(y)   for  the  shock  wave  Is  similarly  elementary  and  an- 
alytic, then  the  equations  (3-9)  m^y  tie  continued  analytically 
Into  the  complex  x-plane  by  Inspection,  l.e,  simply  by  substi- 
tuting for  y  ( =^ )   In  (3«9)  the  complex  variable   I  +  Irj 
That  the  Initial  data  may  be  extended  to  complex  values  Is 
a  crucial  property  that  we  shall  resort  to  later  on  when  we 
solve  the  Cauchy  problem  In  the  complex  domain. 

4 .    Complex  Extension  Applied  to  Equations 
In  Characteristic  Normal  Form 

Thus  far,  we  have  a  formulation  of  the  Inverse  pro- 
blem written  In  two  ways,  namely  as  a  first  order  system  (2.20) 
and  as  a  system  In  characteristic  normal  form  (2.29-31).   In 
Section  3,  we  demonstrated  how  Initial  conditions  behind  the 
shock  are  obtained  from  the  shock  conditions.   The  partial 
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differential  equations  and  Initial  conditions  comprise  the  full 
statement  of  the  Cauchy  problem  for  the  flow. 

We  shall  concern  ourselves  at  present  in  making  sure 
that  such  a  formulation  is  a  proper  one  in  the  sense  that  the 
solution  depends  continuously  on  the  initial  data.   In  the 
supersonic  region,  the  equations  are  of  hyperbolic  type  and 
the  prescribed  initial  conditions  are  sufficient  for  the  exis- 
tence of  a  unique  solution.   However,  in  the  subsonic  region, 
the  equations  are  of  elliptic  type,  which  means  that  initial 
conditions  do  not  specify  a  properly  posed  problem.   Following 
the  procedure  used  by  Garabedian  [2,6,7 ]»  we  shall  take  advantage 
of  the  analyticity  of  the  equations  and  apply  the  method  of 
complex  extension  to  them.   The  resulting  equations  in  the  com- 
plex domain  will  be  of  hyperbolic  type,  so  that  the  Cauchy  pro- 
blem becomes  well-posed  there. 

In  this  section,  we  shall  concentrate  our  attention 
on  the  equations  (2.29-31)  in  characteristic  normal  form.   For 
easy  reference,  we  restate  them  here  as 


^Q  -  ^4-  -\  =  °  '        yp  -  ^_  ^p  =  0  ^ 


(4.1) 


2   2         2   2  /  2   2   2 

T,(c  -u  )u   +  (c  -V  )v   ^  Qc  /u  +v  -c    y   =  0, 


T_(c^-u^)up  +  (c^-v^)Vg  -  Qc  /u^+v^-c^   yg  =  0> 


^p  +  PVXp  -  puyg  =  0   , 


-  24 


where 


Pk       2   2 

_  -uv  ±  c  s/\l    +v  -c 

T+  -       2   2 
c  -u 


We  shall  limit  our  considerations  to  the  subsonic  region  where 

the  Cauchy  problem  is  not  well-posed. 

2   2    2 
In  the  subsonic  region,  u  +v   <  c    and  the  solutions 

of  the  characteristic  equation  (2.23)  are  in  general  complex. 

If  a(x,y)  =  const.   denotes  one  family  of  characteristics,  where 


a(x,y)  =  t(x,y)  +  is(x,y)  , 

then  we   must   have 

p(x,y)    =  t(x,y)    -    is(x,y)    =  const. 

as  the  second  family  of  characteristics  because  the  coefficients 
of  (2.23)  are  real  when  x   and  y  are  real.   Now  a  ~   const, 
and  p  =  const,   m.ay  serve  as  a  suitable  set  of  characteristic 
coordinates  for  the  canonical  system  (4,1)  in  the  subsonic  region 
if  x,y,u,v,^  are  allowed  to  assume  complex  values.   The  system 
(4.1)  will,  however,  be  hyperbolic  only  when  restricted  to  that 
coordinate  system  with  respect  to  which  the  characteristic  curves 
a  =  const.   and  P  =  const.   are  real. 

The  extension  of  the  flow  quantities  to  complex  func- 
tions of  the  complex  variables  x  and  y  may  be  done  by  the 
basic  rules  of  analytic  continuation  because  of  the  analyticity 
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of  these  functions.   Let  us  for  the  moment  keep  x   and  y  real, 
In  the  real  x,y-plane,   t   and   s   are  real-valued  functions. 
Since,  in  addition,  they  are  conjugate  harmonic  functions,  we 
can  perform  the  transformation   (x,y)  — >  (t,s)  .   Then  the 
curves   a  and  P   in  the   x,y-plane  become  the  straight  lines 

a(t,s)  =  t  +  is   ,   P(t,s)  -  t  -  is  , 


in  the   t,s-plane.   For  convenience,  let  us  assume  that  the 
shock  wave  is  mapped  onto  the   s-axis.   We  then  make  a  Judi- 
cious use  of  complex  extension  and  extend  the  variable   s   to 

the  complex  plane  so  that 
s 

^t 

q  =  r  +  is 


-*-  r 


Then  the  corresponding  extensions  of   a   and   P   are 

a(t,q)  =  (t+r)  +  is   , 

P(t,q)  =  (t-r)  -  is   . 

If  we  restrict  our  attention  to  the   t,r-plane,  we  observe  that 
a  and  P   are  real  functions,  i.e.  in  this  plane  the  system 
(4.1)  has  real  characteristics  and  thus,  is  hyperbolic. 
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The  Cauchy  data  for  this  system  Is  now  the  complex 
extension  into  the   r,s-plane  of  the  initial  data  that  was 
defined  for  t   and   s   real.   No  difficulties  arise  from  the 
analytic  continuation  because,  as  we  noted  in  Section  3,  the 
initial  data  for   u,v,p,p   are  computed  from  analytical  formulas 
(3.9  -  10)  which  involve  only  the  basic  operations  of  addition, 
subtraction,  multiplication,  division  and  square  root,  and 
depend  on  a  function  f'(y)   which  is  also  analytic  and  ele- 
mentary. 

The  solutions  we  obtain  shall  be  valid  in  some  section 
of  the   t,r-plane  since   e   is  kept  fixed  as  a  parameter  in 
our  calculations.   If  in  this  surface  of  solution  we  set   r  =  0  , 
then  we  obtain  in  the  t,s-plane  a  line  of  real  solutions  which 
is  of  physical  interest.   To  get  the  solution  in  a  full  region 
of  the  t,s-plane,  we  must  then  solve  the  hyperbolic  problem 
repeatedly  for  a  full  range  of  values  of  the  parameter   s  . 

Since   s   is  a  parameter,  the   t,s-plane  may  still 
be  transformed  conformally  provided  that  the  segment  of  the  s-axis 
which  contains  the  Cauchy  data  is  preserved.   Thus  any  analytic 
change  of  scale  in   s   may  be  introduced.   This  means  that  the 
initial  values  of  any  one  of  the  unknowns,   x,y,u,v,^,   say   y 
can  be  given  by  a  quite  arbitrary  analytic  expression  in   s  . 
Given  a  fixed   s   then,  this  Implies  that   y  ,  as  a  complex- 
valued  function  of   r  ,  can  be  any  complex-valued  analytic 
expression  that  assumes  its  conjugate  value  at   -r  .   In  parti- 
cular  y(r)  may  be  any  continuously  differentiable  function  of 
r  ,  to  which  we  can  approximate  by  analytic  expressions  in   r  . 
This  is  possible  since  we  have  a  well-posed  Cauchy  problem  and 
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any  convergent  limiting  process  on  the  data  still  leaves  us  with 
valid  results.   All  the  other  variables   x,u,v,?//  must  however 
be  initially  related  to  y  by  the  analytic  conditions  of  our 
Cauchy  problem  in  the  original   x,y-plane.   In  essence,  our 
choice  of  the  function  y=y(r)   is  equivalent  to  choosing  a 
path  of  integration  in  the  real   x,y-plane  if  we  were  to  solve 
the  Cauchy  problem  by  successive  approximations.   Since   y(r) 
need  not  be  analytic,  this  path  need  not  be  analytic  either. 
As  we  shall  see,  it  can  have  corners,  cf.  Figs.  2,3. 

So  far,  we  have  indicated  a  way  of  solving  the  Cauchy 
problem  separately  in  the  subsonic  and  supersonic  regions  of  the 
flow.   This  method  is  not  completely  satisfactory  in  the  sub- 
sonic region  because  each  calculation  gives  us  only  a  line  of 
solutions.   Moreover,  in  the  supersonic  region,  this  method 

has  a  serious  defect.   The  system  (4.1)  becomes  singular  along 

2   2     2 

the  sonic  line  since  along  it   u  +v   =  c   ,  i.e.  the  term 


v/u  +v  -c    has  its  branch  point  on  the  sonic  line.   This  implies 
that  characteristics  in  the  supersonic  region  terminate  on  the 
sonic  line  and  form  cusps. 

Let   a(x,y)  =  const.   denote  the  family  of  character- 
istics whose  slope  is   t   and  P(x,y)  =  const.   denote  the 

T" 


second  family  whose  slope  is   x   ,  see  figure. 


a  =  K^ 


shock 


-P  =  const 


body 


sonic  line 
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Let   P-,   be  the  point  where  the  sonic  line  crosses 
the  body.   Consider  the  constant   K  such  that  the  character- 
istic curve   a(x,y)  =  K  crosses  the  body,  any  curve   a(x,y)  =  K  +  e 
for   e  >  0  also  crosses  the  body,  but  the  curve  a(x,y)  =  K  -  e 
terminates  on  the  sonic  line.   If  the  characteristic  curve 
a(x,y)  =  K  hits  the  sonic  line  at   P-,  ,  then  we  know  that  we  can 
solve  for  the  flow  In  the  entire  supersonic  domain. 

However,  if  the  curve   a(x,y)  =  K  terminates  at  a 
point   P   on  the 'body  a  distance  away  from  P,   (see  figure 
below  and  also  Figure  3  in  the  Appendix),  then  we  know  that  a 
section  of  the  supersonic  region  bounded  by  the  sonic  line 
P.,  P   ,  the  body   P-i  P^   and  the  curve   PpP^  which  is  the  limit 
of  the  characteristic   a(x,y)  =  K  +  e  ,  as   e  — >  0  ,  lies  out- 
side the  domain  of  influence  of  our  initial  data  on  the  shock 


wave 


shock 


character- 
istics 


sonic  line 


body 


^11 
It  is  hoped  that  using  the  second  formulation  of  the 

Cauchy  problem  (2.20),  we  shall  be  able  to  solve  for  the  flow 

in  this  hidden  transonic  region.   We  turn  our  attention  to  it 

in  the  next  section. 
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5.    Complex  Extension  Applied  to  a  General  First 
Order  System  In  2  Independent  Variables 

We  shall  concern  ourselves  with  the  system  (2.20)  In 
this  section.   In  matrix  notation,  (2.2G)  has  the  form 


5.1 


Together  with  some  Cauchy  data  along  the  Initial  curve  C,    =  0    , 
we  know  that  we  can  solve  this  system  by  successive  approxi- 
mations wherever  (5.1)  Is  of  hyperbolic  type.   Hov;ever,  as  we 
saw  in  Section  4,  this  is  not  the  case  in  the  subsonic  region. 

By  the  use  of  complex  extension  on  the  independent 
variable   ^  ,  we  shall  demonstrate  how  we  can  obtain  a  system 
equivalent  to  (5.I)  which  is  hyperbolic  in  both  the  subsonic 
and  supersonic  regions.   We  proceed  as  follows:   Consider   ^ 
as  the  real  part  of  the  complex  variable  x       so  that  X    =  ?  +  iri 
Using  the  reflection  principle,  we  can  likewise  extend  all  the 
dependent  quantities   u,v,p,p   as  complex  functions  cf  x.  The 
system  (5.I)  then  becomes 


(5.2) 


U   =  RU   +  S 


Initial  Line 


Note  that  in  general   U^   =  U^  =  -1  U   ,  since  by  the  Cauchy- 
Riemann  equations 


30 


U^  =  (Re  U)^  +  Kirn  U)^  =  (Im  U)^  -  l(Re  U)^  =  -  lU^  , 


Hence  we  may  write   IL.   =  -^  (U^  -  iU  )  .   Similarly,  if  we 

write  X  =   i    -   ±r]    ,    then  U  =  =r   (U*  +  iU  )  =  0  .   We  can  thus 

add  the  term  R*  U_  ■  to  the  right  hand  side  of  (5.2)  obtaining 

X 
the  system 


(5-3)  U^  =  RU  +R*U  +B 


where  the  asterisk  *      denotes  the  conjugate  transpose  of  the 
matrix.   Reverting  to  the  variables   ?  ,  t)  ,  and   i^  ,  equation 
(50)  can  be  written  as 


(5.M  U^  =^^  U^  +^2^^^  +  s 


This  new  system  is  symmetric  hyperbolic  because  the  coefficient 
matrices   (R  +  R*)/2   and   (R  -  R*)/2i   are  Hermltlan  and  hence 
have  real  eigenvalues.   Thus  In  the  new  ^,    x],    t,   -domain,  the 
Cauchy  problem  Is  well  posed. 

Our  derivation  of  the  system  (5.4)  did  not  depend 
in  any  way  on  whether  we  were  In  the  subsonic  or  supersonic 
region.   We  should  then  be  able  to  solve  the  equations  in  the 
transonic  region  without  any  difficulty,   especially  In  that 
section  which  would  have  been  impossible  to  penetrate  by  the 
method  of  characteristics.   The  one  disadvantage  to  this  formu- 
lation is  that  we  have  to  carry  out  our  computations  in  a  three- 
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dimensional  domain  In  order  to  get  a  surface  of  solution  In  the 

real  C>    ^    -   plane. 

From  the  above  discussion,  it  would  seem  as  though  the 
instability  of  the  elliptic  Cauchy  problem,  in  the  sense  of 
Hadamard,  has  been  eliminated  by  the  application  of  a  simple 
complex  extension  to  one  of  the  variables.   This  is  not  the  case, 
however.   We  replaced  the  task  of  solving  an  unstable  elliptic 
Cauchy  problem  in  the  real   C  j  ^  -  plane  with  that  of  having 
to  solve  a  well-posed  hyperbolic  Cauchy  problem  in  the  complex 
C  ,  X  -  plane  plus  having  to  continue  the  initial  conditions  into 
the  complex  domain.   Solving  the  hyperbolic  problem  numerically 
is  a  stable  process.   However,  instability  enters  In  the  analytic 
continuation  of  the  Cauchy  data  into  the  complex  domain  because 
this  process  is  equivalent  to  solving  Laplace's   equation  with 
Cauchy  data.   Hence  in  general,  we  may  not  be  able  to  carry  out 
the  analytic  continuation. 

The  success  of  the  method  in  our  particular  case  de- 
pends on  the  fact  that  we  have  chosen  simple  analytic  curves 
like  hyperbolas  and  parabolas  for  the  shock  shape,  enabling  us 
to  continue  the  Cauchy  data  into  the  complex  domain  in  closed 
form.   Thus  we  readily  obtain  the  analytlcity  of  the  data  in 
the  large  and  we  can  investigate  the  singularities  arising  from 
their  continuation  explicitly. 
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6.    Numerical  Solution  of  the  Cauchy  Problem 

In  this  section  we  shall  present  two  numerical  schemes, 
one  for  the  characteristic  equations  (4.1)  and  another  for  the 
first  order   system  (5.4). 

The  system  (4.1)  Is  of  the  form 


yp,  =  '^_  (x,y,u,v,^)  Xp   , 


(6.1)    ^1  (x,y,u,v,^)  u^  +  B^  (x,y,u,v,V')  v^  =  C^  (x,y,u,v,^)^ 

I    ' 

^2  (x>y.u,v,^)  u   +  B^  (x,y,v,v,i^)  Vg  =  C^  (x,y,  u,  v,  ^/) , 


^p  =  C^  (x,y,u,v,^)   , 


and  we  shall  be  able  to  apply  the  method  of  Massau  [5,13], 
which  is  designed  specifically  for  a  system  of  quaslllnear  hyper- 
bolic equations  In  two   Independent  variables.   The  system 
has  to  be  In  characteristic  normal  form  because  Massau' s  scheme 
takes  advantage  of  the  existence  of  two  distinct  families  of 
characteristics. 

Let   I   be  the  Initial  non-characterlstlc  curve  and 
let   P-,  ,  Pp  ,...,P^   be  a  sequence  of  points  along  I  .   Through 
each  point  on  I  pass  two  characteristic  curves  with  slopes 
1/t   and  1/t_  •   Let   Q  be  one  of  the  two  points  of  Inter- 
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section  of  the  four  characteristic  curves  through  two  adjacent 


points   P.,   and   P^   on   I   such  that   Q  is  on  the  side  wherein 
we  want  to  find  the  solution  x,y,u,v,^'  .   We  can  solve  for 
x(Q)  ,    y(Q)   from  the  finite  difference  equations 

x(Q)  -  x(P^)  =  T^(x(P^),y(P-j_),u(P^),v(P-^),^'(P^))[y(Q)-y(P^)]  , 


(6.2) 


x(Q)  -  ^(P^)  =  ^_[^{V^),y{V^),\x{V^\v{V^),^{?^))[y{Q.)-y[V^)]       . 


Next,  we  can  obtain   u(Q)   and   v(Q)   from 


(6.3) 


A^(P^)[u(Q)-u(P^)]  +  B^(P^)[v(Q)-v(P^)]  -  C^(P^,Q) 


k^{V^)[v.[(^)-u{V^)]    +  B2(P2)[v(Q)-v(P2)]  =C2(P2,Q)  , 


where  x,y,u,v   at   P.,   and   P    plus  the  newly  computed  values 
of  x,y  at   Q  are  used  to  ( 
Finally  for  ^(Q)  ,   we  have 


of  x,y  at   Q  are  used  to  calculate  the  coefficients   A.,B.,C. 


(6.4 


V/(Q)  -  7A(P^)  =  C^(P^,Q) 
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Thus  far,  this  set  of  difference  equations  has  a  trun- 
cation error  of  first  order  because  we  have  effectively  approxi- 
mated the  derivatives  by  forward  differences.   We  can  obtain  an 
error  of  second  "6fder  If  we  use  two  Iterations  to  compute 
x,y,u,v,-^  at   Q  In  the  following  manner.   First  the  equations 
(6.2-4)  are  solved  for  x,y,u,v,^  at   Q,  .   Then  using  the  aver- 
ages |[x(P^)  +x(Q)]  ,  |[x(P2)+x(Q)]  ,  etc.   Instead  of 
x(P  ),  x(Pp),  etc.  In  the  calculation  of  the  coefficients 
T^  ,  A.,B.,C.,  we  use  equations  (6.2-4)  a  second  time  to  obtain 
x,y,u,v,^  at   Q  .   The  combined  effect  of  the  two  Iterations 
Is  equivalent  to  approximating  the  differential  equations  by 
centered  differences,  thus  giving  us  a  truncation  error  of 
second  , order . 

This  scheme  may  be  used  to  solve  the  system  (4.1)  both 
In  the  subsonic  and  supersonic  regions,  the  only  difference 
being  that  In  the  subsonic  region,  all  the  quantities  are  com- 
plex variables  of  r  and  t  .   The  reflection  principle  may 
be  taken  advantage  of,  however,  because   x,y,u,v,''//  are  real 
on  the   t-axls,  and  calculations  need  to  be  carried  out  only 
In  one  quadrant  of  the   r,t-plane. 

To  obtain  a  difference  scheme  with  second  order  accur- 
acy for  the  first  order  system  (5-4),  we  need  to  do  a  little 
more  work.   Let  us  rewrite  i'^A)    as  the  nonlinear  system 

U^  =  A(^,Ti,C,U(^,Ti,C))U^  +  B(e,Ti,C,U(^,r),0)U^ 
(6.5) 

+  c(^,Ti,r,u(?,,ii,r)) 
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with  the  prescribed  Cauchy  data  U(^,ri,0)  .  The  most  obvious 
scheme  to  use  Is  that  Involving  centered  differences.  Super- 
impose a  rectangular  grid  in  the  ^  ,ri,  C-domain  such  that  mesh 
sizes  along  the  C,  11  and  C,  axes  are  h-j^,  h^  and  k,  re- 
spectively.  Then 


U(e,Ti.C+k)  =  U(^,ri,C-k)  +  |-  A(e,Ti,C,U)[U(4+h^,Ti,C)-U(4-h^,Ti,C)] 

(6.6)  +^  B(^,Ti,C,U)[U(e,Ti+h  ,C)  -U(^,T]-h2,C)] 

2 


+  2k  C(^,Ti,C,U)  +  O(h^k)   . 


In  this  form,  we  need  the  values  of   U  at  six  points  to  determine 
U  at  a  new  point.   To  decrease  the  amount  of  computation,  we 
note  that  Instead  of  using  U(^,ri,C)   in  the  computation  of  the 
matrices   A,  B,   and   C,   we  can  use  the  average 

U   =  J  [^{i+h^,-q,0    +   U(4-h^,n,C) 

(6.7) 

I   U(e,T]l-h2,C)  +  U(e,Ti-  hg,^]  -  U(|,Ti,C)  +  0(h2)  . 


2 
This  would  introduce  an  error  of  order   h   in  the  calculation 

of  the  matrices   A,  B,   and  C   because 


A(4,T],C,U  +  0(h^))  =  A(4,Ti,C,U)  +  O(h^)  A„(^,Ti,C,U) 


U' 


=  A(^,Ti,C,U)  +  O(h^)  , 
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If  we  assume  that   | A, J  <  K  in  our  domain  of  solution.   This 

2 
error  of  order  h   in  the  matrices   A,  B,   and   C   does  not 

increase  our  truncation  error  of  order  k   for  U(^,ri,^+k). 

Thus  our  difference  scheme  now  has  the  form 


U(^,ri,C+k)  =  U(^,Ti,C-k) 


+  |-  A(C,Ti,C,U)[U(|+h^,Ti,0  -U(|-h^,7i,C)] 


(6.8) 

+  ^  B(e,Ti,C,U)[U(^,Tl+h2,C)  -U(|,ii-h2,0] 

+  2k  C(^,Ti,C,U)  +  0(k5)   , 

where   U   is  defined  in  (6,7 ).   To  determine   U  at   (t,T],^+k)  , 
we  need  to  know  it  only  at  the  alternate  points   (|+h-,,ri,^)  , 
{ii-h^,r],Cj    ,     (|,Ti+h2,C)  ,    [i,r\-h^,0      and   (e,ri,C-k)  . 

However,  we  note  that  the  values  of  U  at  two  levels 
of  C,      are  needed  in  order  to  compute  it  at  a  third  level.   As 
our  problem  is  stated,  this  is  impossible  since  we  are  given  U 
only  at  the  initial   /^-level.   At  the  start  we  need  an  auxiliary 
scheme  which  will  produce  a  second   il^-level  of   U  values  from 
the  initial   f^-level.   In  order  to  be  compatible  with  our  main 
difference  scheme,  this  a'uxiliary  scheme  must  give  at  least 
second  order  accuracy.   Although  we  do  not  have  to  worry  about 
the  stability  of  the  scheme  since  we  are  to  use  it  only  once, 
we  have  to  be  concerned  with  the  accuracy  of  the  results  on  the 
second   (!;-level.   The  success  of  the  stability  analysis  on  the 
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main  difference  scheme  (6.8),  that  we  shall  make  later,  depends 
on  how  close  to  the  truth  we  are  when  we  assume  that  we  have 
two  levels  of  analytic  initial  data.   Since  the  difference  scheme 
is  an  approximation  to  the  differential  equations,  the  solution 
that  we  shall  obtain  at  the  second  l^-level    differs  from 
the  true  solution  by  some  error.   This  error  will  be  of  an 
oscillatory  nature  since  our  equations  are  of  hyperbolic  type. 
Hence  if  the  mesh  width  in  the   (-direction  is  reduced,  this 
error  function  will  be  damped  out. 

We  shall  now  proceed  with  the  derivation  of  the  aux- 
iliary difference  scheme.   We  expect  beforehand  that  this  scheme 
will  involve  more  than  the  five  points  required  by  (6.8)  in  order 
to  compensate  for  the  lack  of  a  second   (-level.   First  note 
that 


(6.9)  u(e,ii,k)  ^  u(^,Ti,o)  +  ku^:4,Ti,|)  +  o(k^: 


k  2 

This  implies  that  we  need  to  know   Ux^^,!!,^)   up  to   0(k  ) 


Now  from  (6.5) 


U^(^,Ti,|)  ^  A(|,Ti,|,U)  U^(^,Ti,|) 


(6.10; 


+  B(^,Ti,|,U)  U^(^,Ti,|)  +  C(^,Ti,|,U)  ; 


hence  we  see  that  we  actually  need  to  know  U,  U^ ,  U   up  to 


k 


0( k  )   at  the  intermediate  level   ^  =  p  *   Expanding  U,  U^,  U 
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and  using  (6.5)  we  obtain 


U(|,Ti,|)  =  U(e,Ti,0)  +|u^(|,Ti,0)  +  0(k2) 


k 


(6.11)       =  U(e,Ti,0)  +1  [A(|,ri,0,U)U^(e,Ti,0) 


+  B(^,ii,0,U)U^(4,Ti,0)  +  C(e,Ti,0,U)]  +  0(k2)  , 


and 


u^|^,ii,|)  .  u^(^,Ti,o)  +  ^  \J^^^{i,^,o)   +  o(k2) 


(6.12)        +  B(?+h^,ii,0,U)  U(|+h^,ri,0)  +  C(|+h^,ri,0,U) 


-  A(e-h^,Ti,0,U)  U^(|-h^,Ti,0) 


-  B(e-h^,Ti,0,U)  U  (^-h^,7i,0)  -  C{^-h^,r],0,\J)] 


+   O(k^) 


A  similar  expansion  holds  for  U- (.|,ri,p-)  ,      The  expressions  we 
have  obtained  for  U,  U. ,  and  U   involve  only  the  space  de- 
rivatives of  U  on  the  initial  level  and  we  shall  use  centered  dif- 
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.erences  to  approximate  tiiem.   In  computing  the  coefficient 
matrices  A,    B,  and   C   at   (|,  t),  ~  ,    l}[^,-q,:^))      in  equation 
(6,10),  we  shall  use  the  approximation  (6,11)  for  \]{i,r\,:^   )    , 
Altogether,  our  aioxlliary  difference  scheme  Is 

U(|,Ti,k)  =  U(^,Ti)  +  k  {A(?,ri,|,  U)E^ 
(6.13) 

+  B(|,Ti,|,  U)E2  +  C(4,Ti,|,  U)  ]  +  0(k5)  , 


where 


and  where 


U  =  U(|,Ti)  +  P(e,Tl)  , 


+  Th^  [F{i+h-^,r])    »  F(|»h^,ri)] 


^2  =  2n~  [^'(^^^+=^-2^  -  Vi^,^-h^)] 


+  j^   [F{^,r]+h^]    -   F{i,r]-h^)]    , 


F{^,r])    ^  ~-   {A(e,ri,C,U)rTj(^+h.^,n)  »  U(  ^-h-^,!!  )  ]] 

+  2H~  {  B(e,Ti,0,U)[U(|,r,+h2)  -    Vi{^ ,J]-h^)  ]} 
+  C(^,T1,0,U)  . 
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Thus  we  have  a  difference  scheme  which  approximates 
the  differential  equations  to  second  order  accuracy.  Up  to  this 
point  we  have  employed  quite  laborious  means  to  get  a  good 
solution  at   the  second   level.    Next  we  have  to  worry 
about   what    happens  to  our  solution  as  the  calculation  is 
carried  out  from  1^=0     to   C  =  T  .   The  main  point  of  a  stability 
analysis  is  to  examine  the  extent  to  which  any  component  of  an 
initial  function  can  be  amplified  by  Iterations  of  the  numerical 
procedure. 

For  the  purpose  of  this  analysis,  we  consider  A,    B, 
and  C   as  constant  matrices.   We  then  write  the  difference 
scheme  (6.8)  as 


U(e,Ti,C+k)  =  U(e,Ti,C-k)  +  ^  A[U(^-h^,Ti,C)  -U(^-h^,ii,C) 


(6.14) 

-I- 

h 


+  ^  B[U(4,Ti+h2,C)  -  U(^,ri-h2,C)]  +  2kC  . 


In  its  present  form,  we  have  a  three-level  difference  scheme. 
For  convenience  let  us  reduce  it  to  a  two-level  system  In  the 
following  manner: 


U(^,^,C+1^)  =V(e,ri,C)  +  A  ^  [U(?,+h^,Ti,C)  -U(^-h^,Ti,C)] 
(6.15)  +  ^  h"  fU(^,,ri+h2,C)  -  U(?,,Ti-h2,C)]  +  2kC, 


V(^,Ti,C+k)  =  U(?,Ti,a 
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In  essence  we  have  a  system  of  the  form 


(6.16)     ^        {^,r^)    =   R(k)  w'^d,!!)  +  S(k) 


As  our  calculations  go  from  Z    =  0     to   C  =  T  ,  we  are  in  effect 
applying  the  sequence  of  operators   R  ( k)  ,  0  <_  nk  <_  T  ,  on  the 
initial  data   W°(x,y)  ,      As  defined  in  [12],  the  operator   R(k) 


is  said  to  be  stable  if  the  set  of  operators   R  (k)  ,  n  =  l,..=,T/k 
is  uniformly  bounded  for   k   small  enough. 

Following  the  procedure  described  in  [12],  we  assume 

that  (6  =  l6)  is  solvable  for  periodic   W(i;,ri)  ,  which  is  uniquely 

determined  by  the  equations  and  periodicity  conditions.  Then 
we  can  use  the  Fourier  series  representation  of  the  solution  to 
obtain  a  stability  criterion.   We  write   W(f;,ri)   as 

ri(v^+v  T]) 
Z{v^,v^)    e        ^        ^        , 
r-L^r^^l 

where  v  =  2it  ^'-i/Ij-,  and  v  =  2it  ^o/^p  ^^^  "^J^e  frequencies 
corresponding  to  the  periods  L^  and  L„  of  W(|,ri)  in  the 
^-      and      rj-      directions,    respectively,    and      Z(v-,,v    )       is    the 

amplitude    of   the   harmonics    of   the    fundamental    vibration 

27ri(?/L^+Ti/L2) 
e  .      The   ecjuations    (6.1^)    then  become 
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'21    [l^  ^    sm    v^h^    +|-B   sin    v^h^]    1^1 


I  0 

-i(v,x+v_y) 


(6.17)  -  4,c"  ^ 


12 
+   e 


0 


o7   ' 

where   the      G  y:   6     coefficient   matrix 


21[A  T--  £in    v^  h^    +  B  t —   sin   v^h    1  I 

h-,  11  hp  2   2 

(6.18)   G(k,v^,V2^ 

0 


is  called  the  amplification  matrix. 

We  shall  use  the  following  definition  of  matrix  norm: 


(6.19)  Ij  G  II   =  max   i^  =  max   |  Gw 


W  I  ;^0     I  W  I       I  W  I  =1 


/  2  .     ,2 
=  /w 


where   w   =  /w.,  +...+W    .   Now  the  matrix   G  can  be  parti- 
'  '     "  1       m  ^ 

tioned  in  the  form 


G      I 


G 

I 


where   G   is  the   3  x  ^  matrix 


^1  =  2^^^  TT  ^^"^   ^I'^l  +  B  |-  sin  v^h^]   , 


and   1   is  a   ;5  x  5   identity  matrix.   By  construction,  the 
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matrices  A      and  B  are  Hermltian,  hence  their  eigenvalues  are 
all  real.   Let   A„   and  A^   denote  the  maxima  of  the  eigenvalues 
of  A     and  B  .   Then  the  amplification  matrix   G   is  majorized 
by  the  6   y:   6     matrix 


?1  [A^   ^  sin  v^h^  +  ^B  H;  '^"^    "2^2^  ^ 


(6.21) 


and  we  need  only  study  the  eigenvalues  of   p. 

Let   A  ,  A  ,  ...,A/-   be  the  eigenvalues  of  the  matrix 
1  '  c   The  maximum  of   | A . | ,  i  =  1, . . . , 6   is  called  the  spectral 
radius   cr  of   |  '  .   Clearly  <y     jf_  ||  |  '  ||  because  the  maximum 
value  of   I  \~'   w|/|w|   is  not  less  than  the  value  obtained  by 
substituting  for  w  an  eigenvector  of   |  '     corresponding  to  its 
largest  eigenvalue.   Hence  also   (T   ^  II  1  *  II   ^"^^   a  necessary 
condition  for  stability  is  that  there  exists  a  constant   K   such 
that   0-  '^(k,v  ,v  )  <  K  for   0<k<T,   0^nk<^T=Nk. 
Without  loss  of  generality  assume  that   1  <_  K  .   Then  we  have 
that 


^(k,v^,v^)  <k1/>^  -k'^/'^ 


Now  for   k  in  the  interval   0  <  k  <  t  ,  with  small   x 


k/T 
K  ^   <_  1  +  K,  k  ,        K-,   some  constant. 
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hence  a  necessary  condition  for  stability  is  that  the  eigen- 
values of   p  must  satisfy 


A .  I  <_  1  +  0  ( k )    for   0<k<_T  ,  1  =1,...,6 


The  matrix  \^     has  the  triple  pair  of  eigenvalues  of 
the  form 


/ 


(6.22)  A  =  I   (ico  ±  tA-o)^   )   , 


where 


CO  =  2(A^  ^  sin  v^h^  +  Ag^  sin  v^h^ ) 


In  order  that   |a|  =1  ,  co  must  satisfy  the  inequality   |cd|  <  2 
Thus  a  necessary  condition  for  stability  is  that 

(6.23)   |A,  ^  Sin  v^h^  +  A^  ^  Sin  v^h^ |  <  | a,  f-  +  A^  ^  |  <  1 

It  can  be  shown  (see  [12])  that  condition  (6.23)  becomes  suffi- 
cient for  stability  if  the  equality  sign  Is  omitted.   In  that 
case,  the  determinant  of  the  eigenvectors  of   ]  ' 


-  I  |(4-cd2)  -ioovA-oo^  j 


is  bounded  away  from  zero. 

Since  in  our  problem  the  matrices   A   and  B  are  not 
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constant,  but  depend  on  the  solution,  the  quantities   A„   and 
A^  were  computed  at  every  point  of  the  grid  during  the  actual 
calculation  so  that  the  stability  of  tne  run  could  be  checked. 

The  computation  procedures  for  the  coefficients  in 
both  methods  are  essentially  the  same.   Given  u,v ,f     at  a 
point,  the  various  quantities   A(^),  c  ,  p,  p,  A^  {ip)/Aif) 
and   Q  are  calculated  in  the  following  order: 


N^ 


U 


l+f'2(|)       ' 


u  =y7  M   , 


A(^)    =   C^^    [(1-H^)    N?-h2][^l2n2    +    (14^2)^7    ^ 


2       7-1 


2  2 

P        -  ^^ 2 '    JaJJT 


c      =  yA[ilj)    p^  , 


P  =  ^^P  ( TTT    i°e  P^"  ) 


p   =  .A(^)    p^   ^   p 


A'  i>) 


=    -2 


f'(^)f"(|) 

u[i+f'2(i)] 


(1-^l2)n^[^,2n2+(i^2j^_^(^^2)^(^_^2j^2_^2^ 
[(l-H^)N2-tx2][M.V^+(l4vL2)] 


and  finally 


Q  = 


c  A'  [Tp] 

7(7-1)      ATvT"     P 
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To  solve  the  system  (5.4),  we  set  up  a  triangular 
grid  in  the   ^^^l-plane  with  its  base  on  the   ^-axis.   We  are 
able  to  confine  our  calculations  to  the  upper  half-plane  be- 
cause of  the  symmetry  properties  acquired  by  the  functions 
u,v,p,   and   p   under  analytic  continuation.   By  the  reflection 
principle,  any  function   f  which  is  real-valued  on  the  ^-axis 
may  be  analytically  continued  into  the  complex  plane  by  the  rule 
f(^  -  ir|)  =  f  ( ?+irj)  •   We  compute  initial  values  at  every  grid 
point,  then  use  the  auxiliary  scheme  to  obtain  a  second   i!^-level 
of  solutions.   Since  the  latter  is  a  13-point  scheme  using  the 
points  indicated  in  the  diagram,  the  triangular  grid  is  decreased 
by  two  boundary  rows  at  the  second   /^-level.   However,  from  then 
on,  we  use  the  5-point  main  difference  scheme,  which  shrinks  the 
grid  by  only  one  boundary  row  at  each   C-level.   We  have,  in  the  end 
a  tetrahedral  domain  of  answers  in  the   ^,r|,  ^-space,  of  which 
only  its  triangular  section  in  the  real   ^,^-plane  is  of  interest 
to  us  , 
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v.    Numerical  Results 

In  the  preceding  sections  we  derived  a  formulation  of 
the  Inverse  detached  shock  problem  in  terms  of  the  unknowns 
u,v,   and  f    .      We  found  that  the  equations  of  motionj,  treated 
either  as  a  first-order  symmetric  hyperbolic  system  (5.4)  or  as 
a  system  in  characteristic  normal  form  (6.1),  plus  initial  con- 
ditions on  the  shock  constituted  a  well-posed  Cauchy  problem 
in  the  complex  domain.   We  then  devised  difference  schemes  of 
second  order  accuracy  to  solve  the  Cauchy  problem  in  its  two 
versions . 

Our  next  task  is  to  check  the  feasibility  of  the  method 
in  practice.   The  stability  condition  (6.23)  was  derived  under 
the  stringent  assumption  that  the  coefficient  matrices  were 
constant.   In  reality,  our  equations  are  nonlinear  and  we  want 
to  see  whether  the  stability  estimate  we  made  is  still  valid. 
We  also  want  to  know  whether  the  Inverse  problem  can  be  solved 
for  the  flow  quantities  and  the  body  in  a  reasonable  amount  of 
machine  time.   In  addition,  because  of  the  nonlinearlty  of  the 
equations,  there  is  the  question  of  how  much  of  the  body  we  can 
solve  for  in  the  large.   Finally,  although  theory  points  to  the 
three-dimensional  system  (5.4)  as  the  better  of  the  two  versions, 
we  want  to  see  whether  such  a  preference  is  justifiable  in  practice 

All  versions  discussed  in  Section  2  were  programmed  in 
the  FORTRAN  II   language  for  the  IBM  7094  at  the  NYU  AEC  Computing 
Center,   The   u,v,p,p-   and  ^,^^,^, -versions  served  only  as 
checks  on  the   u,v, ^-programs  and  no  significant  results  were 
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obtained  from  theiHo   For  finding  actual  flows,  only  two  programs 
were  used:   Program  I,  which  applied  the  difference  schemes  (6.13) 
and  (6,8)  to  the  system  (5-4)  In  the  three-dimensional   ?,ri^C- 
domaln,  and  Program  II,  which  Integrated  the  system  (6,1)  along 
characteristics  by  Massau's  method  (6.2-4).   Program  I   uses 
8729  machine  locations  and  has  sufficient  storage  space  for  a 
59-lteratlon  run.   For  more  iterations,  tapes  or  discs  would 
have  to  be  used  for  temporary  storage.  A   listing  of  machine 
instructions  for  Program  I  is  included  in  this  paper.   Program 
II  has  two  forms:   the  first  employs  only  real  arithmetic  and 
is  Intended  for  use  in  the  supersonic  region;  the  second  employs 
complex  arithmetic  and  is  used  in  the  subsonic  region.   Program 
Il(real)  takes  up  5524  locations  and  has  sufficient  storage 
space  for  an  897-iteratlon  run.   Program  II (complex)  uses  6718 
locations  and  can  do  a  run  of  at  most  857  iterations.   In  quoting 
the  sizes  of  these  programs,  we  are  including  in  the  count  the 
storage  space  required  by  the  utility  subroutines  used  by  the 
program. 

We  shall  now  discuss  the  results  obtained  for  a  flow 
at  Mach  6  with  the  shock  wave   x  =  3  /1+y   .   In  Figure  1,  we 
exhibit  that  portion  of  the  flow,  body  and  sonic  line  which  a 
single  45-iteration  run  with  Program  I  gave.   This  run  took  I3.2 
minutes  on  the  IBM  7094  and  the  results  were  good  to  four  signi- 
ficant figures.   Table  1  indicates  the  convergence  of  the  answers 
at  a  point  in  the  subsonic  region.   Table  4  lists  the  coordinates 
of  the  interpolated  body  points,  and  the  pressure,  velocities 
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and  Mach  number  at  these  points.   The  ratio  of  the  detachment 
distance   5   to  twice  the  radius  of  curvature   r   of  the  body 
at  the  vertex  was  found  to  be   5/2r   =  .217  .   This  compares 
favorably  with  the  data  obtained  by  Kim  [9],  who  made  experi- 
mental studies  of  flows  past  cylindrical  bodies.  A   second 
quantity  which  Kim  considered  was  the  ratio  of  twice  the  radius 
of  curvature  to  the  sum  of  the  detachment  distance  and  the 
radius  of  curvature.   In  our  case,  this  quantity  was   2r  /&+r 


.82 


In  Figures  2  and  3,  we  exhibit  the  results  obtained 
for  the  same  flow  using  the  method  of  characteristics.   Heavy 
lines  mark  that  portion  of  the  sonic  line  and  the  body  solved  for 
by  Program  II.   In  Figure  3j  we  display  the  section  of  the 
transonic  region  which  lies  outside  the  range  of  influence  of 
the  shock  wave.   Figure  J>   also  Includes  the  continuation  of 
the  sonic  line  in  front  of  the  body  obtained  by  integrating 
the  equations  (6.1)  along  characteristics  in  the  direction  away 
from  the  body.   Finally,  notice  in  Figures  2  and  3  the  three- 
cornered  path  of  answers  obtained  by  using  a  similarly-shaped 
initial  curve  In  the  complex   t,r-plane.    This  illustrates  the 
comment  made  in  Section  4  on  the  possibility  of  using  an  initial 
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path  in  the  complex  domain  which  Is  not  analytic. 

In  the  subsonic  region,  60  steps  were  needed  to 
get  to  the  sonic  line  behind  the  body  with  results  good  to  five 
significant  figures.   Such  a  run  took  1.8  minutes  and  Table  2 
shows  the  convergence  of  the  answers  at  a  typical  point.   Note 
that  to  locate  the  body  and  the  sonic  line  adequately,  four  or 
five  values  of   s   had  to  be  used.   In  the  supersonic  region, 
120  iterations  were  needed  to  get  body  points  with  five  digit 
accuracy.   Such  a  run  took  4,2  minutes  and  Table  3  illustrates 
the  convergence  of  the  answers. 

As  a  check  on  the  accuracy  of  the  calculations,  the 
various  versions  of  the  program  were  used  to  compute  the  flow 
quantities  at  particular  points.   The  answers  obtained  agreed 
with  one  another  to  five  significant  figures  as  Tables  5,  6 
and  7  indicate. 

We  should  like  to  note  that  the  run  exhibited  in  Figure 
1  involves  an  additional  transformation  of  the  x,  (^-domain,  namely 


/.^ 


x'  =  x+mC,   C'^/m  +1   1^,     m  constant. 


This  transformation  rotates  the  integration  path  and  allows  us 
to  approacn  the  body  in  a  direction  which  more  closely  approxi- 
mates the   T    characteristic  directlono   The  best  value  to  use 
for  m  depends  on  the  inclination  of  the   t  -characteristics 
in  the  physical  plane  with  respect  to  the  x-axis.   For  the 
particular  run  in  Figure  1,  we  used  m  =^  1/2  .   This  rotation 
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is  most  needed  In  the  supersonic  region  where  the   T_^-charact- 
erlsticE  slant  toward  the  x-axls. 

Aside  from  the  particular  flow  discussed  above,  we 
considered  various  other  shock  shapes  and  flow  speeds.   As  to 
the  shock  shape   x  =  f(y)  =  f(|)  ,  there  are  some  limitations 
on  our  choice.   Foremost  Is  the  necessity  to  choose  elementary 
and  analytic  expressions  for   f(^)   In  order  that  analytic  con- 
tinuation of  the  Initial  conditions  can  be  performed  In  closed 

form.   A  second  limitation  arises  If  we  consider  the  formulas 

2 
(3.9)  .   The  term  1  +  f  (^)   occurs  In  the  denominator  and 

hence  we  have  to  worry  about  the  zeros  of  this  expression.  For 
example,  consider  the  hyperbola   x  =  /m  -1  /l+y   ,  whose  asymp- 
totes are  the  characteristics  of  the  Incident  flow.   After  analy- 
tic continuation  Is  performed  on   ^  =  y  ,  we  have  that 

f'(:e)  =  v/f?^  [W'iy^^'^  X   , 

where   :«  =  ^  +  Iti  .   To  avoid  the  points   iC  =  +  i/^  at  which 
the  expression   1  +  f '  (x)  =  0  ,  we  have  to  limit  the  extent  of 
our  initial  domain  in  the   ri-dlrection,  which,  in  turn,  con- 
trols the  size  of  the  domain  of  solution  through  the  stability 
condition  (6.25).   The  principal  difficulty  in  practice  is  to 
ensure  that  the  allowable  domain  of  solution  will  not  fall  short 
of  the  body.   Hence  the  choice  of  a  shock  wave   f(jc)   is  also 
Influenced  by  the  locations  of  the  roots cf   1  +  f  (X)  =  0  . 
It  is  desirable  to  have  them  as  sparse  in  number  and  as  far 
from  the  origin  as  possible. 
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While  experimenting  with  different  shock  shapes,  we 
came  upon  a  third  factor  that  has  to  be  taken  into  account  when 
choosing  a  shock.   With  the  use  of  Program  II,  it  was  found 
that  for  M  =  6   the  flow  corresponding  to  the  shock  wave 
X  =  v/m  -1   /1+y   contained  an  envelope  of  characteristics  in 
the  supersonic  region  between  the  shock  and  the  body.   In  terms 
of  the  physical  flow,  this  envelope  could  be  interpreted  as  a 
secondary  shock  occurring  in  the  flow,  a  phenomenon  which  is 
not  unusual.   Although  no  difficulty  was  encountered  in  getting 
to  the  body  and  the  sonic  line  behind  the  body  using  Program  II 
(complex)  in  the  subsonic  region,  we  could  not  reach  the  body 
in  the  supersonic  region  with  Program  II  (real)  because  of  such 
a  singularity  in  the  flow.   Program  I  yielded  even  less  infor- 
mation than  was  obtained  from  Program  II  (real)  because  the 
extent  of  its  region  of  solution  was  limited  by  the  envelope  of 
characteristics  in  the  supersonic  region.   Thus,  we  have  to  con- 
clude that  in  order  for  Program  I  to  give  satisfactory  answers, 
the  flow  behind  the  shock  wave  must  be  continuous  and  free  of 
any  secondary  shocks.   For  the  case  we  considered  at  Mach  6, 


2 

we  were  then  led  to  try  the  hyperbola   x  =  5  \/l+y 

Even  when  Program  I  fails,  however,  we  are  still  able 
to  get  valuable  information  about  the  flow  and  the  body  by  using 
Program  II.   Using  the  complex  arithmetic  version,  the  portion 
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of  the  body  in  the  subsonic  region  can  be  solved  for,  while  the 
real  version  used  in  the  supersonic  region  gives  us  an  indication 
of  singularities  in  the  flow.   For  flows  with  small  Mach  numbers, 
e«g,   M  =  1.2,  we  have  not  yet  found  a  shock  shape  which 
does  not  have  accompanying  secondary  shocks  in  the  supersonic 
region.   Finding  a  shock  shape  suitable  for  Program  I  is  a 
harder  job  for  small  Mach  numbers  because  the  detachment  distance 
5   increases  as  M  decreases.   This  means  that  we  have  to  solve 
the  nonlinear  problem  in  an  even  larger  domain.   Moreover,  from 
the  physical  point  of  view,  secondary  shocks  are  expected  when 
M   is  small. 
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TABLES 


The  following  tables  contain  results  from  the  accuracy 
and  convergence  tests  performed  for  a  flow  at  Mach  6,  7  =  1.4 


/ — 2 
with  the  shock  wave   x  =  3  / 1+y 


Table  1 
Convergence  at  the  point  x   =   3'0696592  , 
y  =  .041458976   using  the  three-dimensional 
method. 


No. 

of 

Iterat 

ions 

^1 

^2 

k 

P 

10 

.05 

.025 

.0125 

43.206604 

20 

.025 

.0125 

.00625 

45.259683 

40 

.0125 

.00625 

.005125 

43.266751 

Table  2 
Convergence  at  a  point  using  the  method  of 
characteristics  in  the  subsonic  region. 


No.  of 

Iterations 

Ar 

X 

y 

p 

15 

.0112 

3.0970604 

.051564741 

35.568517 

30 

.0056 

3.0970272 

.051679900 

35.562432 

60 

.0028 

3.0970185 

.051709075 

35.561677 
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Table  3 
Convergence  at  a  point  using  the  method  of 
characteristics  in  the  supersonic  region. 


No.  of 

Iterations 

Ay 

X 

y 

p 

30 

.008 

3.2109012 

.14448365 

13.024765 

60 

.004 

3-2108428 

.14457851 

13.022444 

120 

.002 

3.2108301 

.14460678 

13.022385 

Table  4 
Interpolated  flow  quantities  on  the  body- 
obtained  from  three-dimensional  scheme. 


X 

y 

P/Pmax 

M 

tan-l  ^ 

V 

3.0657 

0. 

1. 

0. 

0. 

3.0666 

.016902 

.98746 

.14224 

.12363 

3.0687 

.028756 

.96176 

.24144 

.19968 

3.0716 

.040378 

.92476 

.33969 

.28259 

3.0755 

.051729 

.87812 

.^3773 

.35619 

3.0803 

.062780 

.82355 

.53640 

.44625 

3.0860 

.073525 

.76296 

.63603 

.52680 

3.0926 

.083929 

.69838 

.73680 

.60672 

3.1002 

.093971 

.63189 

.83865 

.68614 

3.1088 

.10361 

.56479 

.94329 

.76536 

3.1183 

.11283 

,50107 

1.0462 

.84259 

3.1289 

.12158 

.44406 

1.1438 

-  .91507 

3.1404 

.12988 

.39739 

1.2295 

.97883 
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Table  5 
Results  at  the  point   x  =  3.OO68638  ,  y  =  .024 
after  five  iterations  using  various  three- 
dimensional  formulations  of  the  flow  problem. 


P 

M^ 

u,  v,p,p-version 

42.385768 

.15766106 

i^ff    ,^, -version 

S      Ty 

42.385007 

.15757674 

u,  V , '^-version 

42.584973 

.15757716 

Table  6 
Results  at  the  point   x  =  3. 0591460  ,  y  =  . O72465631  usinj 
two  forms  of  the  u, v , ^-version  in  the  subsonic  region. 


P 

M^ 

No.  of 
Iterations 

3-dimensional  method 
method  of  characteristics 

39.828291 
39.825180 

.25751588 
.25760289 

64 
18 

Table  7 
Results  at  the  point   x  =  3.I388759  ,    ¥   =  .21OO38OI  using 
two  forms  of  the  u, v, ^-version  in  the  supersonic  region. 


P 

m2 

No.  of 
Iterations 

3-dimensional  method 
method  of  characteristics 

21.027724 
21. 022659 

1.6820741 
1,6825560 

56 
36 
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M=  6 


shock        X  =  3  ^\+  v^ 


Figure  I 
Calculation  of  the  flow  based  on  one  run  of  the  three-dinnensionai  method. 


Figure  2 

Portion  of  flow  obtained  using  the  method  of  characteristics.    Characteristics 
are  displayed  in  the  supersonic  region.   The  dotted  curves  in  the  subsonic  re- 
gion indicate  the  intersection  of  the  x,y-plane  with  the  solution  surface  in 
complex  domain. 


shock   X  =  3/1  +y^ 


Figure  3 

Enlargement  of  Figure  2  displaying  the  transonic  region  and  featuring  the 
continuation  of  the  sonic  line  in  front  of  the  shock  wave. 
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FJRMAK/^iX  l<*HGH10    OlMENSiaNS    AHEIJ.2MH    PAINTS     IN    THE    ttS  1 -01  R  EC  t  I  0 
«N/liJX  fl-aN0n.2aM    PBINIS     IN     THE    EIA-DIBtCIlaN// ilXi^HLEFI     END    CBRNE 
*H    0F    L.-<IO    IS   LOCATED   AIF  1?.Ci//i3xHMHl    -  I  PE  1 4 .  T  ,2  >i.HH2    ■  IPEK.r.ax} 
thK     :  IPC  U.r//i5X?HGAMKA     -OPFl*.  I,  I0X2'*H1NI  HAL    H3RUBNTAL    VELBCItV 
K^  lPEll..?//?iX  I6MIMIS    HUN    hiLL     OB  M  .  I  2H     ITERATI0NS.1 

FafiMATI / JiXSTHGIVfJ     SH0CR     HAVE     IS     THE    HYPERHaLA       X     =    S.'SOHTFII     * 
XT>i2)//13X2eH1.4     THIS    CASE     MACh    NLMuER    M    >fK.I) 

F^ftMAM  Ihl/t  )X  lOHSHaCKEH     VII 

F<)fiMan  iHl ) 

FiJHMAII  IOx,  MS.  IP2£20.ai 


,DELXiDELY,OELT.GANHA,Vau.ITEHAI 


KL*3I.BIIJ,KL),D(IJ,KL*5).1J= 


IV  PREPARES     INITIAL     DATA    USING     SH0CK    C0NDIIIk}NS 

SEE    EONS    I  J. 221     IN    PAPER 
SUBRaUT INE    READY 

C0MMaN    N.NX.NY,NT,NX2.NV2, DEL  X.DEL  V.DEL T. GAMMA, H.vau.STi 
C0MM0N    A,B,C.U,U,R,S, I,SUM,HASH.SMACH,I 1 1  ME , SL JPE , GBSH 
CIMENS10N    AU.ai.BI  3,  J»,C(31,0I  J1.HI31  ,S(J}.II  il.SUMIil 
DIMENS  laN    DELXI  II . DEL Y ( II .DEL T 11  I .GAMMA  I  I  I  .M 1 1  I , vaul 1 1 ,1 
DIMENSION    SL0PEI  I  ),MASH< 121,2  I .G0SHI  121.21 
OlMENSiaN    SMACHl 121 ) 
CIMENSIilN    SSI  II.HCAMI  I  I 
ITIME'-: 

GNU- I  GAMMA- (  ).,0.  )  )/IGAMMA*l  l.,O.I ) 
h-l  (Y0U*vau)/(2..O.I  )*(CAMMA/(GA^'MA-I  1.  ,0.  )  I  ) 
00     I     J'I,NY 
NXIIaNX*l-J 

SS(2I-FL0ATF( J-l I'OEIY 
00    I     l-J.NXII 
K-1-NX2 
L-NY*  l-J 
IF(K)     2,2,3 
K-I 
L-J 

SSI  U>IFLaAIFI  I-  l)>DELXI»STAI)Ti 
FPRIME=OER  IwlSSI 
CidS-I  I..0.  1  /SURTFI  I  I. 
SIN-FPH|Me.C0S 
VECIAN=YUU»SIN 
VEC.JBR^Yau«C0S 
tl0ST0N  =  VECNaR*VECNaR 
CEeS0-IGNU«BaST3M*l  I 
VECN0H=C£ESO/VECNa« 

UIK,L  ,  11  =  1 VECN0R*CaSl*(VECTAN<SlN 
UIK.L  ,2  ):<VECT  Arnicas  )-lVECN0H*SIN 
UIK 


x.paw 

61  ,61.3) 


lFPRIMt«FPftlMEM 


*GNLi 


=ss«vau 

I, It.  1 


*l  I  I  l.tO.  )-GNU 


II,)  I-UU.L.3I 


:>  I  I 


U(K 


.  I  I 


£NI 


SMACHl J)=SyflTFlHCAMI I) I 

G3SHI   I,  I  I'UlH.L.  1  I 

GSSHl  1,2  I=UIK.L.2I 

C0Nr INuE 

RETURN 

END 

CaMPuTES    FUNCTION    AI     SECflNU    ^EIi 
DIFFERENCE     SCHEME    ---     SEE     ECN     (i 
SUHR0UT  INE     AIMIIPRINII 
CaMMAN    N.NX  ,NY  ,.'JT,NX2.NV?,DEL  X.DtL 
CBMHaN    A,d,C.U.Q.R.S.T.SUM,HASH,SMj 
CIMENSiaN    All, 31. BIJ, 3). C 111,013), I 
DIMENS  I  arj    DEL  XI  )  I.DELYI  ))  .DEL  TI  I  1  .< 
CI  MEN  SUN    SLaPEI  D.HASHI  12  1,2). G^  Si 
CIMENSIUN    S^-ACHI  121  I 
DIMENSION    Ull 3).U2I 3),IEMP| SI 
CIMEMSiaN    JAH(  )0).GAMI  j) 


GAMMA 

M. 

Yi;u,S 

ME 

iL^PE 

.GJSH 

3) 

II 

1  , 

suMH 

EUUIVALENCE  I JAMI  I  I, JAM!) , IJAMI2) ,JAM2) .IJAM(3) . JAMS), (JAM  1 1 
I.IJA^'l^l.JAMS),  I  JAMI6),  JAM6I  ,  I  JAMiri  ,  JAMTI  ,  |J  AM  I  S  )  ,  JAM8  I  .  I, 
JAMVI  .UAMI   )01,JAM)0).  I  GAM!  21  ,GAM2) 


REHINE 

J 

ITIME= 

1 

NXM):.), 

<-) 

NV)=Ny 

•  1 

Cl=(2. 

0.  liCELX 

C2=7.« 

LtLY 

CdNSI^ 

ELT/l((i.,0.  )*OELX) 

CaNS2= 

[;elt/iii*..o.i 

DELYI 

LJ     7     J 

=I.NY? 

JJ=NXf 

i-J 

KK=J*2 

CD    T     I 

=KK,JJ,2 

CALL     I 

liRLtfi  1  I.J.  JAM 

GAM) 

La     1     K 

-I.N 

Ull«  1  = 

U(  JA-'I.JAMlS.NI-UI  JAM2 

JAM? 

KM /CI 

U2('(  l» 

UI  JA*' J,  JAMrt.tll-Ll  JAMU 

JAM'. 

KH/C2 

U2m*31-(UlJAMl,JAMfl,Ktl)-I&AM2.c 

JAMIi,JAM9,Kt  11)  )  /C2 

CALL     t 

tDIUMI  I, J. SUM 

CALL     I 

tOIUMI 1-2. J, S 

LJ    S    X 

:  I.N 

riK  1=  1 

SUM(K)-SI«(|  1 

CHNSl)*U 

IN) 

CALL     T 

UIIJMI  I-1,J»1 

SUM) 

LALL     I 

DIUMI  I-  1,  J-) 

S) 

ca  <•  ■■ 

I.N 

SUMIK  ) 

I ISUMIK l-SIX 

i*caNS2i 

02  IK 

CALl      A 

'CnARI  1,  J  ) 

CALL     V 

nVEl A.I.U.M.SUM.H \ 

ca   10 

to,  1  )1,  IPHIJ 

MRITE 

APE    '.  I, J,  1  r 

ME  ,  u  ai  1 

.XL  1 

AIIJ.KL 

.  1)  ,d(I J.KLI 

J=I.3I 

KL'  1.1) 

La  ft  K 

I.N 

lEMPIR 

=  U1  JAMS.JAMIO.M'H  JIKI.H(K).CIK)) 

•DELT) 

WRITE 

APE      2, I TtMPlK) . TEMP (K 

3),ll 

I.NI 

IFIJ-I 

7,>),? 

CALL     I 

EIA,  lEMP,  I, J 

0,  21 

CaNHNDt 

REM  INL 

2 

ca  ?   J 

),NV2 

JJ=NXM 

-J 

KK»J*2 

ca  ?    I 

•tK.JJ,2 

in=i-\x2 

JJJ=NY 

-J 

iF(  II  I)  e,tj,i 

111  =  1 

jjj=j 

READ    TAPE    2,1 TEMPIK 

.TEMPIK. 

I.X- 

.N) 

UStS    MAIN    UIFFEHENCE    SCHtMC    SEE    EON    16. dl     IN   PAPER 

SUtlKJl-I  INE     FIRtl  IPRINI  I 

CUMMaN    1,.^X,NY,Nl.NX2.NY2.0£LXtDElY,OELT.GAMM.H,vau.STA>tTX,pai 


CaMMaN    A,B,C,U,U.K, S( T, SUM.NA 
1  CIMENSiaN    AI3,3).HI 1,3I.C( 3). 

I  CIMENSUN    DElxI  D.DELYI  l),DEL 

I  CIMENSiaN    SiaPEI  I  l,MaSMM2l.2 

CIMENSiaN    S^'ACHI  121  J 
i  CIMENSiaN    TtMPI )l 

CIMENS  taN    J  INI  )?l,l>I,'4l  3) 

EQUIVALENCE    IJINin.JIND.UI 
X).IJlNISI,JINbl.(JI'4l6),JlN6) 
XJ1N9).IJINI  I0I,J1NI0I.1GINI2) 
I  CI=D£LI/neLX 

1  C2=ueLT/LELY 

1  C3-12.,0.1»C£LI 

CB     t     KLM=),NI 

1TIME=KLM. 1 
Ny2KL»'  =  NY2-»<LM 
ca     3    J=  I  .NY2KLM 
JKLM= JiKLM 
JJ  =  NX-  l-JKLM 
KR«?»JKLM 
La     3     l=KK,JJ.2 
CALL     b^RDtHt I. J. JIN.GINI 
ca    I    kH.N 
I  CIK  l  =  LI  JINI.JIN6.K  )-Lil  J1N2.JIN 

RIK|ibIJIN3,JlNd,Kl-UUINb,JlN 
I  flIR»3)=UIJINJ,JItja,K*il-IljlJ2» 

CALL     At^CI  1.  J,  ).  ),  1,-1,  JIN, ^ni 
CALL     MQVEI A.Q,S.U,H,U) 

ca  ra    16,7 ).  iPHiM 

6     MRItE  lAPE  .M,  J,  IIIME  .(  (A  I  IJ, 
XJ='I.3),KL=I.31 

r         ca  2  RO,N 

I        lEMPlK  1:UI J  INS, J  IN  10. K  1 

12  UIJIN5,J1NIC,K)>IEMP|K) 

If (J-  n    3, 5, 3 
•b  CALL    Tt-EIA,  TEMP,  I,  J,:, 2) 

3  CONTINUE 

U  CALL    lL.IPUr 

HEM  INI      1 
M£TUH\ 
END 

C 

CABCUAR      CaMPutES  MATRICES  A.a.C 
SUHR0LT1NE  ABCBARI I  I, 121 
CBMMBh  N.NX,NY,Nt,NX2,NY?,l)EL 
CBMMBh  A,B,C,U.U.R.S,  I,SIJM,HA 

1       CIMENSiaN  AI3. 3),B1  3. 3),CI  1) . 

I  OMENS  laN    OELXI  I  l.OELYl  )  )  . Ut I 

I  ClMtNSiaN    SLBPEI  1),mASMM2I.2 

CIMENSiaN    SMACHl I?)) 
CIMtNSUN    KaBKI  IO).ijaaK|  3) 

EguivALENCE   iKBBKi  n  .KaaK  n.i 
xui.Kaii'KK),  UBBKi  s),KaBKb),iKa 

X .KBBKt  )i IKaaKI 91 .KSJKV) , IKa^K 
I  C0NST  I<I2.. C. )*DEL  X 

CaNST2^2.«D£LY 
1  CaNSI  .•'DELI/12., a.  ) 

CALL   eaHL'ERi  1 1,  i2,KaaK,G0aKi 

C0     I    K3  I.N 


SH,SMACH,I  TIML.SL.: 
CISl.KISl .St3l ,  l<  3 

I  .OAMHAI  1  I   ,H  M  ) 

J Sril 121.21 


M2I  .JIN<!)  ,  I  JINI  \  I 

ij  I N 1  n  .  J 14  n  ,  ( J I 

GIN2I 


I  JI  Id,  JIN'^.K*!) 


•  SlNl*L?»mK)»C3«C  IX 


i2.i)Eir/2i   FaH  Sue 


x.n 

L Y.DELI ,G 

MMA.H. 

Y^U.STA 

Sh. 

i-ACH.lIi" 

.SLJPt 

,(.aSH 

u(  1 

,K|31,SI3 

.IM1  . 

SUMlil , 

II  1 

.^AMMAI I ) 

HI)  1  ,Y 

i^l\)  ,l> 

l.i.JSMI  12), 21 

KaaK(2i .KaaK? 

.  i».a^i^ 

1 31 .K0a 

aKI6).Ka0K6) , 

ptaaKi/ 

1  .KaaKf 

(  10 

,KaaKiO). 

k;jaKi2 

I  ,^aaH2 

PAGE  b 

CIKI'lkJIKBaK  l,K0^K6.K  l-Ulr(00K2.KS0K7,KI  I  /C  BNSI  I 
SIKIi  lU(K0aKt,Kl30Kd.K  l-U(KadKli,Ke0KQ,Kn/C0NST2 

CALL     ArtC(  I  I  .  12.  I.  !•   ■•  I  .Kt)0K,&00KI 
CALL     HAtfEI  A,O.R,e.S.QI 


ca  2  < 


:  I.N 


-CiNSTl. 
CALL     Add  1  1. 
RETURN 
tND 


:hu 


*UlKMC<Ktl* 


12.  I.  1.1.2,1 


CI6CIUM  cuts    NECtSSARY    BUT 

C  C0MPUIES      A«U-SUB-i 

C  STBRES    IT     IN    ALL 

SUBaaUIlNE     TEOJUMl  UJ,«LL  I 

C0MM(JN    ^,N)l,Ny,Nl,NK2.NY2,D£LK,0tLy.UtLT,GAMN»,H.Yau,START<,i 

CaMMBN    a,B,C,U.Q.R.S.I.SuM,»ASH.S"ACH,l  IlHE,SLi!P£.&0SH 

I  tlMENSLJN    aU.3I.BI3,3l,CI3l,Ulil.«IJ>.SIil.riil.SUM(3),U(6I, 


I  .21 


ILI  91 


t  CIME.-SldN    OELKl  I  I.UELYI  I 

I  ClHLNSIiJN    SL0PEI  I  ).wAShl 

CIMENSIJN    S"ftCHI  121 ) 

1  CIMENiliJN     ALLl  31 

LlHE^SIkJN    L  I  T01.FLI  31 

EUUWALENCE     IL(1).Ln.lL 

X  I, LA)  ,iLi  r  ).Lr),  iLidi.m 

X.fLJl 
C)-7.*CELX 
C2-2.»CCLV 
CALL    bCKCtm  1*1.J.L,FL  I 

Ca       I      K: I.N 

C;iKI'lulLl.L6.K)-UIL2.L'.K))/CI 
ljim5l  =  lFLI*ILHLl.L6,K.i|-u<L2.L'.l'* 
HIKI-lu(Li,L8,KI-ulLU,L'*,K>WC2 
1  HU»JI^HFLi«U(L5,LB.K*3ll-IEL2»l.lLU 

CALL     fl«Cl  1*1, J.  1.1.  I.O.L.fll 
CALL     "mUEM.U.S.B.R.UI 

L^  ft  K  =  l,^ 


H(  I 


VdUI 1) .P0MI I > 


ALL(K  I  =  S1K  I 

«ETUflN 

END 


•Q(K 


*Cln  1 


C^fl'UTES    THE    HERHITIAN    MATRICES       A    * 

B    =     1«    -    «     STAH)/2I        AND    C 
SEt    EQN    (5.1*1     IN    PAPER 

AJO    C    AT    (L  1>L2) 
A     IS    N^I    C*)MPUieO    IF    l3 
B     IS    ,*3T    CaHPUTED     IF    LU 
MPUIEO     IF    LS 


I 


IS   C0MPUTED    IF  Li  ■■ 

IS   CBNPUIED    IF  LU  =    I 

IS   CflMPUIED    IF  L5  =    I 
SuBrtJLf  INE    ABCfL  I,L2.L3.LU,l5.L6.maC«,TRUCK1 

CUMMiN     ,^,NX.\y,NI,NK2.NT2,DElX,0£LY,OELr.&A»MA,H,yau,STAKTK. 

Zi^fMU    A.B.C.U.U.B.S.l  ,SUM.haSh,SMACm,IT1m£  .SL0PE,L.JSH 

LIM£NSIJ\    AH,J(,BI3.tl.CI5l,&l3I.RI3I.SI3l,ri3I.SUMI31.UI61 

L  l»E.<SIJN    OELXl  n.OELYI  n  .DELTM  UOAMMAI  I  )  .H  1 1  I  .YBUI  11  ,P0WI  I 

CIME^Sl^N    SL0PEI  ll,>tASHII2  1.2I.GidSHI121.21 

CIME.JSIJN    SWACMl  121  I 

CIHENSIJN     ItMPI  lJ,aC0NrKl3,3) .blGAI 3,3) tXPRlll  ,  YPft I  I  1 , TPRI 1 ) 

CIMtrtSlJN    HACK  I  10),  TRUCK!  i),NlX|  IO),FIXt  31 

EQUIVALENCE    INixl1),NlX)HNIX(21.NI«2l,(N|xi5l,NIxil.(NIx(U 


Xl,lNIXlSI.NU5).INIXl6I.NIX6I.INlX17l,Nix;i,lNlKI61.NIXB),lNlX(9 
XNIX9).lNIXIIO),NlxlOI,lFIKI1).fIx1l,IF|xi2l,flX2l,IFIx(31,FU3l 

C0    IS     1-1,10 

NIXI  I  )>maCk( 1 1 

C0    16     1*1.3 

FIXI I I^IRUCKI II 

ZZ    9    K=  I.N 

IF(L6  I      l,r,2 

TEMPI K  |:UIN1X5.NIX10>K  t 

IEHPIKt3l-F  tXUU(NIXS.NlXlO,K*9) 

C0  Ta   'i 

lEMPIK  )  =  R(K  ) 

GO  la   V 

TEMPI  K  |:IUIMX  UNIX 

TEMPI  K.  3  1-1  IFIXUIUINIXl,NU6.t*5>*U{NIx2,NIx7,K*Sm»IFIX3»UINlx3 
«,NIXB,xt3ll»lfIX2»U(NlXl.,NU9,K«JI)l/U. 
CaNT INU 


»U(NIX2,N1X7.K)*UIN|X3,N[X 


*U(MX(i,NlX9t 


IFIL6- 
mR=1 


B.e. 13 


G0    10     U 
13         INK =3 
U         CALL    IFEIBIGA, IEMP.L1,L2.L5.INR) 

C0    10     1=1. N 

Ca  10  J'I,N 
1       AC0NIRI 1,JI>BIGAIJ,I ) 
n    AC0NTRI  I,J»31»-AC0NIHI I,J*3) 

|F(L)  )i,  3,i« 
It     C0  11  l'l,N 

Ca  11  J=1,N 

111  AII.J)=IBIGA(l,J)*ACaNTRII,JI)/(2..0.) 
3  IFILUlS.S.ft 

6  C0    12     I'l.N 

CB     12     J^I.N 

112  B(  l,J  )-IBIGAI  I,J  I-AC0NTR1  I ,J1  I  /lO.  .2.  I 
S  RETURN 

END 
C 
CB0HOER  APPLIES    REFLECT10N    PRINCIPLE    F0R    B0UNOAHT    GRID    P0INTS 

suBKaur  iNt   BaROER  im  i,M2.N00K,H0aii  i 

CaMMjJN    N,NX,NY,NT,NX2.NY2,0£L  X  .DEL  Y  .OELT  .GAMMA  ,H  ,  Y  £u  .  S  I  AR  I  X  . 

C  IMENS  UN     Li;0K  I  6  I  .ItailK  I  31  ,N00KI  10) 

NTUNY.l 


La»K(  1 1 

L00)(|  !\ 


KI  ?  |, 


Ml* 


L0tfKI U 

L00KI Sl  =  M2-  I 

LaaKl 61=M2t  1 

ca     I     I\G=0,( 

IFUauKI  INCH    2.2, 

i.e0Ki  1 4Gi'2-La0Ki 

t-BBKl  I^G-JI=-l.0 

C0    10     I 

hBaKl  lNG-3)-I.O 

CBNTINUt 

N0BKI  I  )  =  LeaKl31 

NBBKI 2)=LB0KI2) 

N00K(  31>L00K( I) 


NBBKI hl<LB0KI I ) 
\BBKI SI'LBBKI  I) 
\BaKI  61=Lk)(?KI4) 
NBkltl  n-LB0Klitl 
NBBKI  eULKBKIfil 
N00KI 9  I^LBBKI $1 

N0BK1  i:  i=LBeK<<«i 


C2    I 


1  iG  = 


Kl  IN&1-'JX2 
|FMIN,:Bk1    Ii,1i,S 

NBBKI   |i4G)=M  IMBBK 

NflUKl  I  ,Lt5l=Ny  l-NBiJKl  ING* 

CBNI INUt 

RETURN 

END 


t  C/MPuTtS    THE    SPtCTl 

SIBKES    II    IN   OESl S 

C>-ARACIER    BF     The    I 
SUBKDIIINE    CtASClfl.*D, DESIST 

CIMtNSUf*     ANOI   i,  3) 
BEE  =  ANCI  1.  1  l*A'gOI2,21 

Aft(,^i  i'-EE«eeei-i'i.»iiANoi  i. 
.5)) 


RADIUS  JF  THE  3X1  i 
TAKING  AOVANTAGF  Bl- 
HATRIX   AND 


Al^DI  1  ,21  (AND!  1.2  1 


IFI ARL  1 
CESISt= 
C»  TB  1 
WBBI'Si. 


1.2.2 
-11111. 


HIFI AKGI 

»SF{  Ib66*«i;JT 
ftBBr2:  JrtSFl  IBEE-RiJBt 
RBBTi=aBSF[ flNOI J 
CES  ISI ^"H"  IFI RiZ 
Rkrtj4(. 

tNO 


t)BI2,'<BdI3l 


MAlw IX    VECTBK    «ULT IPLIER 
SUBRBliI  INt     MAtfeU1,/2.fi./l*./'..^6  1 
ClMt-JSUN    ;!lJ.3),;2l3l,/JMl,/(*M,il.;S(ll,/^(il 
CALL    "  WIEI  71,Z2.MI 
CALL    ynvitHit,lb,l(>i 
RETURN 
END 


E  SCAPECBAT    FUR    WAV 

A'^D    STZR6S    PfirfUUCT     IN    VE 
SUSRBUIINE    "AVIEIUI.vP.Ul) 
CBMMaN    N,NX,.MY,Nl,NX2.Ny?.0tL 
CIMENSIBN    U  II  3,  3l,a2I  31  ,l^iM) 
Li    /    J^ 1, 


riPLICS    I 


,UtLY,UELI,GAMMi,l 


ALk 


MC.O, 


D.C) 
J.I  )•' 


Ur  SPLhS    BUT     THE    ANS-ERS 

SUbRBLllNfc     auTPuT 
CBMMBN    N,NX,NY,NT,'*X2,NY2.0£LX.DELY,UELT,GAMMA,H,yBU,STARlx 


PAGE  a 

CBMMON    A.B.C,U.0,R.S,T.SUM,NAS^.SMACH,II1HE.SLBPE 
ClMENSIBNAIl, I),B(3.i).CI3),Ut3),K(3),SI3l,Tl3l, 
CIMENSIBN   OELXl  ll.DELYI  1)  .OELH  1  )  ,o4MMA|  I  )  , 
C  IMENS UN    SLBPEl  1  l.xASHl 12  1,2 t ,GBSHI 12  1 ,21 


Jul  I  I .PBWI  I 


CIMENSl. 


aCmI  1211 


BUNCH 

kRiTE     ^UIPUI    TAPE     6.3 

IFI  II  l-E-1  I     I..S.S 

INCRE= 1 

11=  I 

JJ  =  NX 

G0     10     ( 

INCRE=; 

11=  IT  l''E»2 

JJ  =  NX-  I-!T1NE 

CB     I     1  =  1  I,  JJ.  INCRE 

YY^l I  FLBAIFI  I-l  )«DEL« 

xx  =  auNC^••Fu^crl  YYi 

WRITE     JUTPUI     lAPt     6,2 

XACHI  11,1  GiiSH  I,J*2l.  J 
RETURN 

FBRMATl/l 3X,  1P7E  la. 7/ 
FBRMAIl  ////  lUX  IHM,   \Ji 

XHMACH     r,UM8EH//  I 
END 


•OELTj/SURTFI t SLBPE'SLBPEI 


lARIXI-ISLBPE* 

.y Y.  lOBSHi 1 .J) 


H    AND    STBRES    I  I 
R 
KACH    NL1M8ER     AT      I 


CBfPuIES    CBEFFlCltM    VAlRI 
SEE    EUN    12.2C1     IN    PAP 

BNLY    CBMPUieS    PRESSURE     A^D 
MFEN   L  7    =    2 

ALSJ    CBMPUTES    VECIBH    S    A'^U    SIBRES    11    IN   C    IF    L'l 
SUBRBUT  INE    THE[E.rEMP,Ll,L2,LS.L7) 

CBMMBN  N,NX,'JY,;jI,NX2,NY2.DEL  x,l)tLy,DtLI  ,&AMMA,.|,Y 
CHMMBN  A,B,C,U,Q,«,S,T,SUM,«ASh.S''4CH.  ITIMfc  .SUPF. 
CIMENSUN  A|3,i),B[5,3),CI3l.l<l31.-<l31,St3),rMI,S 
C  IMENS  ur*  OELXl  1  l.DELYI  M  ,0£L  H  1  1  ,^AMHA(  I  I  ,H(1  I  .YJ 
CIMtrjsUN  SLBPEl  I  I.hASHI  12  1,2  I  ,ij^SHl  I2).<;i 
CIMENSUN    SMACMl  121  I 

Ll>*tNSl,;N    £15,  3  I,  TEMPI  31  ,mCAM1  1)  ,1)11  I  ,C2U?I  1  I  ,L1  (  I 
CIMENSIBN    SlalEl  1I,S0UND1  1  I  ,  uSli  I  I  I  ,  VSU  I  II  .Ci  (  I  I 
BUGGER  =  S(.RIF  I  I  SLBPE"  SL  LlPt  )  •  I  1 .  ,0.  )  ) 
GB    10    I  2e,2'i,i01,L  7 
PL4NT  =  FlOAIH  lIlMt  1 
CB    TB     '.\ 

PLANT  =FlBAIf (  I  I  IME-  I  I 
GB    IB     _M 
PLANI=(0.5,C.0I 

SLANI=1SLBPE»PLANT»D£L II/ULGUER 
C=TEMP( J1/Y2U 
FP«AYM:CtR|VlD) 
FPPHyM:DL'£R  IVIO) 

Clll=llFiaAIFILl-ll«OeLXl»SIAHIXI-SLANr 
CI2l  =  FiaAIFIL2-l  l»DELY 
FPRIME^DEH IVIOI 
Cl  =  l  l.,O.I«lFPRAYM.FP«AYM| 
IFICH  111     1^.  16.  15 
IKCll  ()  )     IS  17,  l'l 
WKirt     JUTPUt     TaP£     6,  18,1  l,L2,rPKAYMl  I  I  .F-'HAYMI?) 

CB    IB    li 


HPWtSSUML ,9X1 


^}    fZ".   ^'uTPur 

IS    NB.7trfC 

^u.Sia"*'  «,Pk;w 

U-'ISI  .Ul^l  ,61,  31 
Ol  I  >  ,>'JMl  1  ) 


PA&E 


to 


.  TEMPI  6)  .F  P'tAV*'!  I  )  .F(>ltAYMl2) 


CZ-IGAWMA-I  I..0.  I 
CJ=I  I  .,0.  UC? 
C>«=(  1  -.0.  1-C? 

FuN=( yju«YaLi/ci 

C5'(FljN»C'.l-C2 
IFICSI  III     r,  I  I, 7 
1FIC5I 21)    r,  13, 7 
WRITE    OUTPUT    TAPE 
G0    I^    ^3 
C6=(C2«FUN)»C3 
OSO^IfPI  ll'IEMPI  I  I 
VSO'TEt-PI  JI.TEMPI2I 

tSU  =  OSl.*VStl 
SIATE=£XPFILiCFIC6/FLN)»i 

IFI  STAIE(  1)  I  W.20,  19 

IFI STATEI2)  I  19.21.  19 

MHITE  dUTPUI  TA^t  6.22.L1,L2 

a  T^   23 

t^n^<iM"=^  (H-(Qso/i2.,o.  1  n•lCA^»'A-^  I.  ,0. )  1 1  /UAKKA»siftTE) 

S0UNO-Caf"'a*$TATE*RHi]GAM 

1FIS0UNO1  n  I    2U,2S.2ii 

IF(  SidijNOI  2)  I    2>t.26.2x 

MRITE    kJUlPUI    TAPE    6,27.L  l,L2.LSa(  I  )  rUSUI2l  .VSO(  I  I  .VSQI2) 

C0    Ta 


b.  1W.L  1.L2, TEMP 


AMNAKC5 


RH0  =  EXPF(L0CFIRHt)GAI 

02    T0     I  12, 1C. 121  ,L  7 

WASM( L  l,2l=STATE»RHaGAM.l 

>.ASH(  L  li  n>IEMPI  3) 

FCAH=tSO/SaLNO 

iHACHlL  1  1=SGRIFIHCAN(  1 | I 

CeSHl I  I.  II^IEMPI  I  ) 

CdSMC  L  l.21  =  tCH»'l  21 

GJ    10     c 


IGAM: 


(1.,0.l 


C2v?=SdUNC-VS0 

C2U2=5.JUNC-US0 

YUVI  =  T£Mp(  1  i.itupi 21 

El  1.21 

=-C3v2 

I6RM=- 

FPRIME.El t.2l 

E(2. >  1 

iC2u2 

C-C2U2 

*H2..0.  (•YUVl'FPRlMbUI 

FPRIHE* 

lER"! 

EU.il 

=-SLJP£/BU&GE« 

1:12,21 

:TI:(<M-ISLaPE»D) 

El  1.  1  ) 

=EI2,2)*I12.,0. 

IFIDI  1 

I)    2,  1,2 

IFIDt  2 

))    2,<<.2 

WHITE 

OUTPUT     TAPE     &.S 

ITEMPII  1 

.TEMPll 

•31,1=1 .3 

CALL     E 

«IT 

CUUG'C 

•OUCGER 

L0     1     1 

=  1.2 

z^    1   J 

-1,2 

ti i.j  ) 

^el  I, JI/OBUG 

ce  6   1 

=  1,2 

tl3.1  1 

MCC.I 

EII.il 

:|  0..C.  1 

IF1L5I 

a.e.9 

HULK=- 

I 12.,C.1»FPRATM 

FPPftYM. 

(CU«FUN 

•C6I-(Gam 

C5-C6I 

tUEUE= 

1  S2UN[:»HULK*RHi9 

/f&AHMA- 

(GAMMA- 

1.  .0.  1  1  1 

CM  1=-! 

UUEu£»(yUVI»TEBM) )/06UG 

Ci21-  lUUtUE»UYUVl«FPRlMt1»C2L21  1/D9U& 

C131--ITe*'Pl2t««Mii;i/BUCGEK 

AETURt. 

F0RMAII // I0II9HOEIERMINANT    IS    7E<ie/ IOX6E20.A/1  0X2  lb) 

FBRMAll// 10<ICMC5    IS    iE«0/ IO»2IS,i.E20.ii1 

F0RMAI 1// 101 lOHC 1     IS    /ERB/ t0l21S,2E20.8l 

FiaRHAT(//IOxUHAlPSI)     IS    itH0/lO)>2IS) 

FiR'*ATI//lOKl2HC»-2    IS    ZER0/215,UE2O.B1 

END 


FUNCT  IJH    CERIVIARGI 

CERIV-H3..C.I»ARG1/SURIF(I  I  .  ,C.  1  •  1  AHG»*RG  I 
RETURM 
END 


CDQERIV 

1  FUNCT  Ufi    ODERIVI  ARGI 

I  C0O-1  l.,0.)*tARG*A(lGI 

1  CENeH-SgRTF(C0CI 

I  CDCRIV»li.,C. l»OENflM»( 1 1..0. 

I  RETURN 

END 


-I IARG*ARG)/C0G1 I 


CFUNCT 


FUNCT  IBN  FUNCT(ARG) 
FUNCT-?. -SaRrFI  l.*IARG*ARGn 
RETURN 
END 


